阅读说明：本技术 一种基于双向积分的卫星重力场反演方法和系统 (Satellite gravity field inversion method and system based on bidirectional integral ) 是由 肖云 王丽兵 成伟 赵玉 于 2020-05-25 设计创作，主要内容包括：本发明涉及一种基于“双向积分”的卫星重力场反演方法和系统。其中,基于“双向积分”的卫星重力场反演方法包括：获取待反演卫星重力场中轨道的初始历元和末尾历元；获取卫星重力场反演模型；将所述初始历元和所述末尾历元输入至所述卫星重力场反演模型完成所述卫星重力场的反演。本发明提供的基于“双向积分”的卫星重力场反演方法和系统,通过采用卫星重力场反演模型对卫星重力场进行反演,能够大大减弱系统误差对卫星重力场反演结果的影响,提高卫星重力场反演的精度。(The invention relates to a satellite gravity field inversion method and a satellite gravity field inversion system based on bidirectional integration. The satellite gravity field inversion method based on the bidirectional integral comprises the following steps: acquiring an initial epoch and a tail epoch of a track in a gravity field of a satellite to be inverted; acquiring a satellite gravity field inversion model; and inputting the initial epoch and the tail epoch into the satellite gravity field inversion model to finish the inversion of the satellite gravity field. According to the satellite gravity field inversion method and system based on the bidirectional integral, the satellite gravity field inversion model is adopted to invert the satellite gravity field, so that the influence of system errors on the satellite gravity field inversion result can be greatly weakened, and the satellite gravity field inversion precision is improved.)

1. A satellite gravity field inversion method based on bidirectional integration is characterized by comprising the following steps:

acquiring an initial epoch and a tail epoch of a track in a gravity field of a satellite to be inverted;

acquiring a satellite gravity field inversion model;

and inputting the initial epoch and the tail epoch into the satellite gravity field inversion model to finish the inversion of the satellite gravity field.

2. The method for satellite gravitational field inversion based on "two-way integration" of claim 1, wherein the process of constructing the satellite gravitational field inversion model specifically comprises:

acquiring measured data of a gravity satellite; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data and inter-star range data; the orbit data includes: an initial epoch and a last epoch of the track;

correcting the initial value of the forward orbit of the orbit by using the initial epoch as an initial integral value by adopting a forward numerical integration method to obtain a corrected initial value of the forward orbit;

correcting the initial value of the reverse orbit of the orbit by using the tail epoch as an initial integral value by adopting a reverse numerical integration method to obtain a corrected initial value of the reverse orbit;

taking the corrected initial value of the forward orbit as an integral initial value, and obtaining a forward reference orbit by adopting a forward numerical integration method;

determining a forward reference inter-satellite distance variability according to the forward reference orbit;

taking the corrected initial value of the reverse orbit as an initial value of integration, and obtaining a reverse reference orbit by adopting a reverse numerical integration method;

determining a reverse reference inter-satellite distance variability according to the reverse reference orbit;

and constructing a satellite gravitational field inversion model according to the forward reference orbit, the forward reference inter-satellite distance variability, the backward reference orbit and the backward reference inter-satellite distance variability.

3. The method for inverting the satellite gravitational field based on the bidirectional integration according to claim 2, wherein the constructing the satellite gravitational field inversion model according to the forward reference orbit, the forward reference inter-satellite distance variability, the backward reference orbit, and the backward reference inter-satellite distance variability specifically comprises:

constructing a forward observation equation according to the forward reference orbit and the forward reference inter-satellite distance variability;

determining a forward satellite gravity field inversion model according to the forward observation equation;

constructing a reverse observation equation according to the reverse reference orbit and the reverse reference inter-satellite distance variability;

determining a reverse satellite gravity field inversion model according to the reverse observation equation;

and determining a satellite gravity field inversion model according to the forward satellite gravity field inversion model and the reverse satellite gravity field inversion model.

4. The method for inverting the satellite gravity field based on the bidirectional integration according to claim 1, wherein before the obtaining of the initial epoch and the last epoch of the orbit in the satellite gravity field to be inverted, the method further comprises:

collecting actual measurement data of a gravity satellite in the gravity field of the satellite to be inverted; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data and inter-star range data; the orbit data includes: an initial epoch and a last epoch of the track;

and preprocessing the collected actually measured data of the gravity satellite.

5. A satellite gravity field inversion system based on 'bidirectional integration', which is characterized by comprising:

the data acquisition module is used for acquiring an initial epoch and a tail epoch of a track in a gravity field of a satellite to be inverted;

the model acquisition module is used for acquiring a satellite gravity field inversion model;

and the gravity field inversion module is used for inputting the initial epoch and the tail epoch into the satellite gravity field inversion model to finish the inversion of the satellite gravity field.

6. The system of claim 5, further comprising:

the measured data acquisition module is used for acquiring measured data of the gravity satellite; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data and inter-star range data; the orbit data includes: an initial epoch and a last epoch of the track;

the forward orbit initial value correction module is used for correcting the forward orbit initial value of the orbit by using the initial epoch as an integral initial value and adopting a forward numerical integration method to obtain a corrected forward orbit initial value;

the reverse orbit initial value correction module is used for correcting the reverse orbit initial value of the orbit by using the tail epoch as an integral initial value and adopting a reverse numerical integration method to obtain a corrected reverse orbit initial value;

the forward reference orbit determination module is used for taking the corrected forward orbit initial value as an integral initial value and obtaining a forward reference orbit by adopting a forward numerical integration method;

the forward reference inter-satellite distance variability determining module is used for determining forward reference inter-satellite distance variability according to the forward reference orbit;

the reverse reference orbit determination module is used for obtaining a reverse reference orbit by taking the corrected reverse orbit initial value as an integral initial value and adopting a reverse numerical integration method;

the backward reference inter-satellite distance variability determining module is used for determining backward reference inter-satellite distance variability according to the backward reference orbit;

and the satellite gravity field inversion model building module is used for building a satellite gravity field inversion model according to the forward reference orbit, the forward reference inter-satellite distance variability, the backward reference orbit and the backward reference inter-satellite distance variability.

7. The system of claim 6, wherein the module for constructing the satellite gravitational field inversion model specifically comprises:

the forward observation equation building unit is used for building a forward observation equation according to the forward reference orbit and the forward reference inter-satellite distance variability;

the forward satellite gravity field inversion model determining unit is used for determining a forward satellite gravity field inversion model according to the forward observation equation;

the reverse observation equation building unit is used for building a reverse observation equation according to the reverse reference orbit and the reverse reference inter-satellite distance variability;

the reverse satellite gravity field inversion model determining unit is used for determining a reverse satellite gravity field inversion model according to the reverse observation equation;

and the satellite gravity field inversion model determining unit is used for determining a satellite gravity field inversion model according to the forward satellite gravity field inversion model and the reverse satellite gravity field inversion model.

8. The system of claim 5, further comprising:

the data acquisition module is used for acquiring the gravity satellite actual measurement data in the gravity field of the satellite to be inverted; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data and inter-star range data; the orbit data includes: an initial epoch and a last epoch of the track; and the preprocessing module is used for preprocessing the acquired actually measured data of the gravity satellite.

Technical Field

The invention relates to the technical field of satellite gravity field inversion, in particular to a satellite gravity field inversion method and system based on bidirectional integration.

Background

The earth gravitational field is the basic physical field reflecting the distribution, motion and changes of earth materials. The method can accurately determine the fine structure of the earth gravitational field and the space-time change thereof, and has extremely important functions and significance in the research of military national defense construction, national economy, space science, earth science and related subjects. The study of the earth's gravitational field is a fundamental task and is also the core and hot spot of research in the field of geodety. Satellite gravity measurement opens up a new era of human exploration of the earth's gravitational field and is considered to be a highly efficient gravity exploration technique with the greatest value and application prospect in the beginning of the 21 st century.

The earth gravity field is inverted by utilizing gravity satellite data, namely a satellite gravity field inversion model is constructed, and the earth gravity field inversion model is a core task of an earth gravity satellite. Common methods for inverting the earth gravity field using gravity satellite data are classified into time domain methods and space domain methods. The time domain method mainly comprises the following steps: the method is suitable for solving CHAMP (ChallengMinisatellite Payload) satellites and GRACE (Gravity Recovery and Climate Experiment) satellites. The airspace method is suitable for resolving the earth gravity field from GOCE satellite (gravity field and step-state Ocean Circulation Exploore, gravity field and static Ocean current exploration satellite) data. Among the many satellite gravitational field solution methods, the kinetic method is a widely adopted classical method.

When a dynamic method is adopted to solve an earth satellite gravity field inversion model, the solution of a Newton motion equation and a variational equation is involved, the orbit integration is required to be firstly carried out, then a method equation is constructed based on residual information, and the gravity field is further solved. The numerical integration is to directly perform numerical integration on the satellite motion equation and the variation equation, and gradually obtain the satellite position, the satellite speed and the state transition matrix at any time by taking the satellite position and the satellite speed of the reference epoch as initial values.

The numerical integration method gives an approximate value on a discrete step point, and inevitable errors, namely a truncation error and a rounding error besides an initial value error exist. In either the single-step method or the multi-step method, errors (including initial value errors and rounding errors) generated in a certain step propagate, the errors are accumulated, and the corresponding numerical method is stable only when the accumulation of the errors is controlled.

Orbit integration is actually the process of extrapolating future satellite orbits by integrating the equations of satellite motion based on a set of initial values. During the extrapolation process, the longer the estimation interval, the larger the error accumulated by the integration, due to the influence of error propagation. The one-way error accumulation effect of the orbit integration makes the normal equation information obtained based on the error propagation law worse, resulting in larger deviation of the estimated bit coefficient.

Therefore, aiming at the situation of one-way transmission of system errors, the technical problem to be solved in the field is to provide a satellite gravity field inversion method which can greatly weaken the influence of the system errors on the satellite gravity field inversion result and improve the inversion precision of the gravity field.

Disclosure of Invention

The invention aims to provide a satellite gravity field inversion method and a system based on bidirectional integration, so as to greatly weaken the influence of system errors on a satellite gravity field inversion result and improve the satellite gravity field inversion precision.

In order to achieve the purpose, the invention provides the following scheme:

a satellite gravity field inversion method based on bidirectional integration comprises the following steps:

acquiring an initial epoch and a tail epoch of a track in a gravity field of a satellite to be inverted;

acquiring a satellite gravity field inversion model;

and inputting the initial epoch and the tail epoch into the satellite gravity field inversion model to finish the inversion of the satellite gravity field.

Optionally, the process of constructing the satellite gravity field inversion model specifically includes:

acquiring measured data of a gravity satellite; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data and inter-star range data; the orbit data includes: an initial epoch and a last epoch of the track;

correcting the initial value of the forward orbit of the orbit by using the initial epoch as an initial integral value by adopting a forward numerical integration method to obtain a corrected initial value of the forward orbit;

correcting the initial value of the reverse orbit of the orbit by using the tail epoch as an initial integral value by adopting a reverse numerical integration method to obtain a corrected initial value of the reverse orbit;

taking the corrected initial value of the forward orbit as an integral initial value, and obtaining a forward reference orbit by adopting a forward numerical integration method;

determining a forward reference inter-satellite distance variability according to the forward reference orbit;

taking the corrected initial value of the reverse orbit as an initial value of integration, and obtaining a reverse reference orbit by adopting a reverse numerical integration method;

determining a reverse reference inter-satellite distance variability according to the reverse reference orbit;

and constructing a satellite gravitational field inversion model according to the forward reference orbit, the forward reference inter-satellite distance variability, the backward reference orbit and the backward reference inter-satellite distance variability.

Optionally, the constructing a satellite gravity field inversion model according to the forward reference orbit, the forward reference inter-satellite distance variability, the backward reference orbit, and the backward reference inter-satellite distance variability specifically includes:

constructing a forward observation equation according to the forward reference orbit and the forward reference inter-satellite distance variability;

determining a forward satellite gravity field inversion model according to the forward observation equation;

constructing a reverse observation equation according to the reverse reference orbit and the reverse reference inter-satellite distance variability;

determining a reverse satellite gravity field inversion model according to the reverse observation equation;

and determining a satellite gravity field inversion model according to the forward satellite gravity field inversion model and the reverse satellite gravity field inversion model.

Optionally, before the obtaining of the initial epoch and the last epoch of the orbit in the gravity field of the satellite to be inverted, the method further includes:

collecting gravity satellite actual measurement data in the gravity field of the satellite to be inverted; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data and inter-star range data; the orbit data includes: an initial epoch and a last epoch of the track;

and preprocessing the collected actually measured data of the gravity satellite.

A satellite gravity field inversion system based on 'bidirectional integration', comprising:

the data acquisition module is used for acquiring an initial epoch and a tail epoch of a track in a gravity field of a satellite to be inverted;

the model acquisition module is used for acquiring a satellite gravity field inversion model;

and the gravity field inversion module is used for inputting the initial epoch and the tail epoch into the satellite gravity field inversion model to finish the inversion of the satellite gravity field.

Optionally, the system further includes:

the measured data acquisition module is used for acquiring measured data of the gravity satellite; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data and inter-star range data; the orbit data includes: an initial epoch and a last epoch of the track;

the forward orbit initial value correction module is used for correcting the forward orbit initial value of the orbit by using the initial epoch as an integral initial value and adopting a forward numerical integration method to obtain a corrected forward orbit initial value;

the reverse orbit initial value correction module is used for correcting the reverse orbit initial value of the orbit by using the tail epoch as an integral initial value and adopting a reverse numerical integration method to obtain a corrected reverse orbit initial value;

the forward reference orbit determination module is used for taking the corrected forward orbit initial value as an integral initial value and obtaining a forward reference orbit by adopting a forward numerical integration method;

the forward reference inter-satellite distance variability determining module is used for determining forward reference inter-satellite distance variability according to the forward reference orbit;

the reverse reference orbit determination module is used for obtaining a reverse reference orbit by taking the corrected reverse orbit initial value as an integral initial value and adopting a reverse numerical integration method;

the backward reference inter-satellite distance variability determining module is used for determining backward reference inter-satellite distance variability according to the backward reference orbit;

and the satellite gravity field inversion model building module is used for building a satellite gravity field inversion model according to the forward reference orbit, the forward reference inter-satellite distance variability, the backward reference orbit and the backward reference inter-satellite distance variability.

Optionally, the satellite gravity field inversion model building module specifically includes:

the forward observation equation building unit is used for building a forward observation equation according to the forward reference orbit and the forward reference inter-satellite distance variability;

the forward satellite gravity field inversion model determining unit is used for determining a forward satellite gravity field inversion model according to the forward observation equation;

the reverse observation equation building unit is used for building a reverse observation equation according to the reverse reference orbit and the reverse reference inter-satellite distance variability;

the reverse satellite gravity field inversion model determining unit is used for determining a reverse satellite gravity field inversion model according to the reverse observation equation;

and the satellite gravity field inversion model determining unit is used for determining a satellite gravity field inversion model according to the forward satellite gravity field inversion model and the reverse satellite gravity field inversion model.

Optionally, the system further includes:

the data acquisition module is used for acquiring the gravity satellite actual measurement data in the gravity field of the satellite to be inverted; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data and inter-star range data; the orbit data includes: an initial epoch and a last epoch of the track;

and the preprocessing module is used for preprocessing the acquired actually measured data of the gravity satellite.

According to the specific embodiment provided by the invention, the invention discloses the following technical effects:

according to the satellite gravity field inversion method and system based on the bidirectional integral, after the initial epoch and the tail epoch of the orbit in the satellite gravity field to be inverted are obtained, the inversion of the satellite gravity field can be completed by adopting the satellite gravity field inversion model, so that the influence of system errors on the inversion result of the satellite gravity field can be greatly weakened, and the inversion precision of the satellite gravity field is improved.

Drawings

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings needed to be used in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings without inventive exercise.

Fig. 1 is a flowchart of a method for satellite gravitational field inversion based on "two-way integration" according to an embodiment of the present invention;

FIG. 2 is a schematic diagram of the concept of "two-way integration" in an embodiment of the present invention;

FIG. 3 is a flowchart of a gravity field inversion method for GRACE for low-orbit gravity satellites according to an embodiment of the present invention;

4 a-4 c are track residual diagrams before correcting the initial data value of a GRACE track of a 24H long arc segment according to an embodiment of the present invention;

5 a-5 c are corrected track residuals after "bi-directional integration" of 24H long arc GRACE track data in an embodiment of the present invention;

6 a-6 c are track residual error graphs before correcting the initial data value of the GRACE track of the 6H short arc segment in the embodiment of the present invention;

7 a-7 c are corrected track residuals after "bi-directional integration" of 6H short arc GRACE track data in accordance with embodiments of the present invention;

fig. 8 is a schematic structural diagram of a satellite gravitational field inversion system based on "two-way integration" according to an embodiment of the present invention.

Detailed Description

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

The invention aims to provide a satellite gravity field inversion method and a system based on bidirectional integration, so as to greatly weaken the influence of system errors on a satellite gravity field inversion result and improve the satellite gravity field inversion precision.

The invention provides a satellite gravity field inversion method and system based on bidirectional integration, which has the following basic ideas: when the orbit numerical integration is carried out, the inverse integration is added, and the two integrals with opposite directions are utilized to more accurately integrate the orbit, so that the inversion calculation is carried out to obtain higher precision.

In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in further detail below.

Fig. 1 is a flowchart of a method for inverting a satellite gravity field based on "two-way integration" according to an embodiment of the present invention, and as shown in fig. 1, a method for inverting a satellite gravity field based on "two-way integration" includes:

and S1, acquiring an initial epoch and an end epoch of the orbit in the gravity field of the satellite to be inverted.

And S2, acquiring a satellite gravity field inversion model.

And S3, inputting the initial epoch and the last epoch into the satellite gravity field inversion model to complete the inversion of the satellite gravity field.

The construction process of the satellite gravity field inversion model specifically comprises the following steps:

and acquiring the actually measured data of the gravity satellite. The measured data includes: an initial epoch and a last epoch of a track.

And correcting the initial value of the forward orbit of the orbit by using the initial epoch as an initial integral value and adopting a forward numerical integration method to obtain a corrected initial value of the forward orbit.

And correcting the initial value of the reverse orbit of the orbit by using the tail epoch as an initial integration value and adopting a reverse numerical integration method to obtain a corrected initial value of the reverse orbit.

And taking the corrected initial value of the forward orbit as an integral initial value, and obtaining the forward reference orbit by adopting a forward numerical integration method.

And determining a forward reference inter-satellite distance variability according to the forward reference orbit.

And taking the corrected initial value of the reverse orbit as an initial value of integration, and obtaining the reverse reference orbit by adopting a reverse numerical integration method.

And determining a backward reference inter-satellite distance variability according to the backward reference orbit.

And constructing a satellite gravitational field inversion model according to the forward reference orbit, the forward reference inter-satellite distance variability, the backward reference orbit and the backward reference inter-satellite distance variability.

The construction process of the satellite gravity field inversion model mainly adopts a construction idea based on bidirectional integration. The principle of "two-way integration" is schematically illustrated in fig. 2, in which forward integration is performed by taking the initial epoch of the track as an initial value, time is moved forward by one step, integration is performed once, and the correction is performed according to the initial condition. The inverse integration takes the track end epoch as an initial value, reverses the time backward by one step and corrects the time according to the initial condition. In the "bidirectional integration" example in fig. 2, in order to demonstrate the effect, the present invention intercepts an arc segment of an orbit for 100 minutes as a real orbit, and performs forward and reverse integration by taking two ends of the arc segment as initial epochs of forward and reverse integration, with an integration step of 5 minutes and an integration time of 50 minutes. Due to the accumulation characteristic of the integration error, the effectiveness of the integration step and the integration interval is a main factor influencing the integration precision, and practice shows that the integration step of 5 seconds has better integration precision. Therefore, in the subsequent gravity field inversion process, the method preferably selects 5 seconds as the integration step length.

In order to further improve the accuracy of the gravity field inversion, the process of constructing the satellite gravity field inversion model according to the forward reference orbit, the forward reference inter-satellite distance variability, the backward reference orbit, and the backward reference inter-satellite distance variability specifically includes:

and constructing a forward observation equation according to the forward reference orbit and the forward reference inter-satellite distance variability.

And determining a forward satellite gravity field inversion model according to the forward observation equation.

And constructing a reverse observation equation according to the reverse reference orbit and the reverse reference inter-satellite distance variability.

And determining a reverse satellite gravity field inversion model according to the reverse observation equation.

And carrying out weighted fusion processing on the forward satellite gravity field inversion model and the reverse satellite gravity field inversion model to obtain a satellite gravity field inversion model.

The following describes a specific process of constructing a satellite gravity field inversion model by using a KSG integrator as an example.

When the KSG integrator is used for bidirectional numerical integration, the positive and negative integration formulas have difference in sign.

The whole process of constructing the satellite gravity field inversion model is as follows:

in the satellite gravity field inversion solution, the parameters to be estimated include: the dual-star state vector (position, velocity), the dual-star accelerometer bias and scale factor, and the gravity field model bit coefficient. Recording the inversion resolving order of gravity field as N, and assuming the position vector of the satellite asVelocity vector of

Is a vector to be estimated, an

Wherein the content of the first and second substances,

x, y and z respectively represent three axes of a coordinate system and are three components of the position. v. of_x,v_y,v_zIs a velocity three component. b_x,b_y,b_zThe deviation is three components. k is a radical of_x,k_y,k_zIs a scale three component. C_nm、S_nmAll gravity field model bit coefficients, n is the model order, m is the model number, when m is 0, coefficient C_n0Called the band harmonic coefficient, when n ═ m, coefficient C_nm、S_nmCalled the fan harmonic coefficient, when n ≠ m, the coefficient C_nm、S_nmReferred to as the field harmonic coefficients.

Then the corresponding state vectorInitial value of vector to be estimatedPartial derivative phi of, andinitial value of vector to be estimatedPartial derivatives of

Comprises the following steps:

assuming the integration step size is h, productThe divider order is i (a 14-order KSG integrator is actually used, i is 14). Because the KSG integrator starts by adopting a central iteration initialization process (which is not the key point of the patent and is not specifically described), the integral starting needs to calculate i node valuesWherein

The integral vector isAndrespectively carrying out forward integration and reverse integration, wherein the calculation formula of each node in the integral starting stage is as follows:

1) when in forward integration, selecting the initial epoch of the known observation arc section as the integral initial epoch, and recording the initial time as t₀The initial value of the integral vector is recorded as

Andthen:

wherein the content of the first and second substances,the right function is integrated for the start time,is t_kThe vector of the integral of the time of day,

is t_kIntegral right function of time (i.e. t)_kPerturbed acceleration experienced by the satellite at time of day).

2) When reverse integration is carried out, the last epoch of the observation arc section is selected as the integration starting time, and the initial time is recorded as t₀The initial value of the integral vector is recorded asAnd

then:

it can be seen from the above derivation that the coefficients in the KSG formulas of the forward integration and the backward integration are completely the same, and h in the forward integration formula can be replaced by-h in the implementation process, i.e. the backward integration formula.

The partial derivative of the initial value of the parameter to be estimated of the reference orbit and the instantaneous state vector can be obtained through the numerical integration, the orbit residual error and the inter-satellite distance variability residual error can be calculated and obtained by utilizing the satellite orbit data and the inter-satellite distance observation data, and the observation equation of the parameter to be estimated of the single epoch (comprising the initial state parameter, the accelerometer deviation and the scale parameter and the gravity field position coefficient) can be established through the following formula after the orbit residual error and the inter-satellite distance variability residual error are obtained:

Δr_iis the (i ═ A or B) orbital residual of the satellite i,Is the inter-satellite distance-variability residual error,

is an inter-satellite distance-variability reference value, r_i0,B_i0,K_i0Respectively, initial state parameters, accelerometer bias and scale parameters for satellite i (i ═ a or B), β are gravity field coefficients,solved in the process of integrating the orbit motion equation,

the partial derivative of the inter-satellite distance variability to the initial state parameter and the gravity field coefficient is determined byAnd (4) linear combination.

The two observation equations can be arranged to obtain respective error equations, and the single epoch error equation is uniformly expressed in the following form:

wherein, V_iFor observation errors, L_iFor orbital residuals or inter-satellite range-rate residuals, A_iThe error equation coefficient matrix is formed by combining the partial derivatives, the delta X is the variation of the parameter to be estimated, and the value of i is from 1 to the number of single arc section epochs. In practical calculation, a plurality of epochs are usually used as a set to form a matrix and the matrix form of the error equation of the multi-epoch is as follows:

V＝A·ΔX-L (12)

where V is the error vector, L is the residual vector, and A is the coefficient matrix. The equation can be derived from the error equation as:

wherein A is^TIs the transpose of the coefficient matrix A, N is the left matrix of the equation, B is the right matrix of the equation, and W is the weight matrix of the parameter to be estimated (usually taken as the identity matrix).

In actual calculation, according to the characteristics of parameters to be estimated, a normal equation is expressed in blocks:

x_L、x_Grespectively representing local parameters (initial state parameters of A star and B star, accelerometer deviation parameters) and global parameters (accelerometer scale parameters of A star and B star, gravity field coefficients), N₁₁，N₂₂Are sub-matrices corresponding to local parameters and global parameters,is the right vector of the normal equation.

After the orbit observed value GNV1B data of A star and B star are processed by the error equation (formula 12), the orbit observation equation (formula 9) can form a normal equation

Similarly, each arc intersatellite speed measurement observation KBRR (K-band range-rate) can form a normal equation

In combination with the orbit data and the KBRR data, there is the following relationship:

wherein the content of the first and second substances,

σ_kbrrfor KBRR observation accuracy, σ_orbAnd the track observation precision is obtained. In the calculation process, the observation precision of the KBRR is 0.2nm/s, and the observation precision of the track is 2 cm.

The normal solution is adopted for solving the normal equation, the local parameters of the normal equation are eliminated, and the normal equation about the global parameters can be obtained as

The above formula is denoted as

Solving the equation of the reduced law to obtain the global parameter variable quantity of

Will be Δ x_GBy substituting the following formula, the local parameter variation quantity Deltax that can be solved_L

The above is a forward solution for a single arc segment. In practice, according to the observation characteristics of the satellite, in order to obtain a high-precision solution result, a plurality of arc data are generally accumulated and solved, a time-varying gravity field is generally solved in a month unit, and a static gravity field needs to use observation data of a longer arc.

In order to ensure the accuracy and the integrity of the data, the method further comprises the following steps before the acquisition of the initial epoch and the last epoch of the orbit in the gravity field of the satellite to be inverted provided by the invention:

collecting gravity satellite actual measurement data in the gravity field of the satellite to be inverted; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data, and inter-star range data.

And preprocessing the collected gravity satellite measured data, specifically preprocessing the measured data (the measured data comprises the initial epoch, the tail epoch, accelerometer observation data, star sensor observation data, inter-satellite range observation data and the like) such as gross error elimination, down-sampling and the like, and completing data preparation.

As another embodiment of the present invention, a low-orbit gravity satellite GRACE is taken as an example to specifically describe the satellite gravity field inversion method based on "two-way integration" provided by the present invention.

The inversion process of the gravity field of the low-orbit gravity satellite GRACE by adopting the satellite gravity field inversion method based on the bidirectional integration specifically comprises the following steps:

step 1: four main load 1B-level data (mainly orbit data, inter-satellite distance measurement data, accelerometer data and star sensor data) of a low-orbit Gravity satellite GRACE (Gravity Recovery and clearance Experiment) are collected. The acquired four main load 1B-level data are preprocessed (including gross error elimination, error correction, space-time standard unification, interpolation and vacancy filling, down sampling and the like), and two sets of data sets of 24H and 6H (namely, arc lengths are 24 hours and 6 hours respectively) arc sections are formed. Each set of data contains 1B-level data for the four classes of principal loads, orbit data, inter-satellite range data, accelerometer data, and star sensor data, GNV1B, KBR1B, ACC1B, and SCA1B, respectively. The two sets of data represent different arc lengths for comparing the two-way integration effect at different arc lengths. The track error accumulates differently for different arc lengths. Theoretically, the longer the integrated arc length, the larger the error accumulates, and accumulates quickly.

Step 2: and (3) taking the initial epoch of the GNV1B data of the GRACE satellite orbit arc segment obtained in the step (1) as an initial value, carrying out forward integration by arc segments by adopting a KSG integrator (5 s is taken as a step length, and the default integration step length is 5s without special description later), so as to obtain a forward reference orbit and a partial derivative, and carrying out initial value correction of the orbit arc segment initial epoch by arc segments.

The orbit initial value correction only uses an orbit observation equation, and the parameter to be estimated is only a satellite state vector, and the specific correction process is as follows: 1) the track residual is obtained using the forward reference track and the track. 2) Only the partial derivatives with respect to the state vector are extracted from the partial derivatives, forming a coefficient matrix. 3) And constructing an orbit law equation by using the orbit residual error and the coefficient matrix. 4) And solving a law equation to obtain the correction quantity of the initial value of the track. 5) And correcting the initial track value by using the initial track value correction quantity.

And step 3: and (2) taking the last epoch of the GNV1B data of the GRACE satellite orbit arc segment obtained in the step (1) as an initial value, carrying out backward integration by arc segments by adopting a KSG integrator, and carrying out initial value correction on the last epoch of the orbit arc segment by arc segments, wherein the correction process is consistent with the correction of the initial value of the forward orbit, and only the observation data is inconsistent (the observation data are all GNV1B data, but the initial epoch of the GNV1B is taken as an initial value for integration when the forward initial value is corrected, and the last epoch of the GNV1B is taken as an initial value for integration when the backward initial value is corrected).

And 4, step 4: and performing forward integration by taking the initial epoch of the corrected arc section as an initial value to obtain a forward reference orbit and a partial derivative, and constructing a forward reference inter-satellite distance variability based on the forward reference orbit and the corresponding orbit, wherein the construction process is to use the dual-satellite forward reference orbit for difference to form the reference inter-satellite distance variability.

And 5: and taking the tail epoch of the corrected arc section as an initial value, performing reverse integration to obtain a reverse reference orbit and a partial derivative, and constructing a reverse reference inter-satellite distance variability based on the reverse reference orbit and the corresponding orbit, wherein the construction process is consistent with the construction process of the forward reference inter-satellite distance variability.

Step 6: in the fusion resolving stage, two fusion resolving strategies can be adopted, and the method specifically comprises the following steps:

a) bit coefficient fusion solution

And constructing a normal equation by using the forward reference orbit and the forward inter-satellite distance variability, and acquiring a forward satellite gravity field inversion model according to the normal equation. And constructing a method equation by using the reverse reference orbit and the reverse reference inter-satellite distance variability and calculating to obtain a reverse satellite gravitational field inversion model. And performing weighted fusion on the two groups of models to obtain a fused satellite gravity field inversion model. The whole bit coefficient fusion solving process is shown in fig. 3.

The weighted fusion strategy adopts a spectral combination method based on the gravity field model bit coefficient. The gravity field position coefficient is estimated by adopting two modes, and a linear unbiased estimation model of the gravity field position coefficient is expressed as follows:

wherein the content of the first and second substances,

for post-fusion model bit coefficients, β₁For the forward estimated bit coefficients, β₂For the backward estimated bit coefficients, W₁Estimating the weight of the bit coefficient for the forward direction, W₂Is the weight of the bit coefficients estimated backwards. According to the least square principle, W is taken in consideration of the fact that the same set of observation data is used in two ways₁＝W₂＝0.5。

b) Method equation fusion solution

And (3) constructing a normal equation by using the forward integral orbit and the forward inter-satellite distance variability, constructing a normal equation by using the reverse integral orbit and the reverse inter-satellite distance variability, performing weighted fusion on the two sets of normal equations, and resolving to obtain a fused satellite gravity field inversion model, as shown by dotted lines in FIG. 3. The model refers to a set of data which is actually a data file and contains the gravitational field spherical harmonic coefficient C_nm、S_nmIn the actual application process, the model data file is used as input, and different business processes can be carried out.

In order to verify the specific effect of the satellite gravity field inversion method based on the bidirectional integration, the bidirectional integration inversion calculation experiment is carried out by adopting GRACE satellite data in 12 months in 2005. The numerical integrator adopts a 14-order KSG integrator (KSG is a fixed-order, fixed-step, linear and multi-step integrator suitable for a second-order differential equation, the names of which are obtained from the last names of three scholars researching the numerical integrator: Krogh, Shampine and Gorden), the integration step length is 5s, the integration arc length comprises two groups of 24H and 6H, and a mechanical model adopted in the orbital integration process is shown in the following table 1:

TABLE 1

Test results

The forward integration and the backward integration are performed on the GRACE track data of the 24H long arc segment, the track residual before the initial value correction is shown in fig. 4 a-4 c, and the track residual after the iterative correction of the initial value of the track is shown in fig. 5 a-5 c. Forward integration and backward integration are performed on the GRACE track data in the short arc segment of 6H, the track residual before initial value correction is shown in fig. 6 a-6 c, and the track residual after iterative correction of the track initial value is shown in fig. 7 a-7 c. The above results show that: the initial value correction can effectively reduce the initial orbit error, the forward integration and the reverse integration can both obtain a reference orbit with higher precision, and the one-way accumulation of the system error in the integration process is effectively reduced, so that the transfer of the error in the gravity field inversion process is reduced, and the gravity field inversion precision is improved.

Compared with the prior art, the invention has the following advantages:

1. and the bidirectional integral is provided, so that the accumulation of the track integral unidirectional error can be effectively avoided, and the initial track resolving precision is improved.

2. The bidirectional integral method is applied to the gravity field resolving, so that the accumulation and diffusion of integral errors can be effectively inhibited, the system errors are eliminated, and the inversion precision of the gravity field is improved.

In addition, for the above-mentioned method for inverting the satellite gravity field based on "two-way integration", the present invention also provides a system for inverting the satellite gravity field based on "two-way integration", as shown in fig. 8, the system includes: the system comprises a data acquisition module 1, a model acquisition module 2 and a gravity field inversion module 3.

The data acquisition module 1 is used for acquiring an initial epoch and a last epoch of a orbit in a gravity field of a satellite to be inverted. The model obtaining module 2 is used for obtaining a satellite gravity field inversion model. The gravity field inversion module 3 is configured to input the initial epoch and the last epoch to the satellite gravity field inversion model to complete inversion of the satellite gravity field.

In order to improve inversion and accuracy, the system further comprises: the device comprises an actual measurement data acquisition module, a forward orbit initial value correction module, a reverse orbit initial value correction module, a forward reference orbit determination module, a forward reference inter-satellite distance variability determination module, a reverse reference orbit determination module, a reverse reference inter-satellite distance variability determination module and a satellite gravity field inversion model construction module.

The measured data acquisition module is used for acquiring measured data of the gravity satellite. The measured data includes: an initial epoch and a last epoch of a track. And the forward orbit initial value correction module is used for correcting the forward orbit initial value of the orbit by using the initial epoch as an integral initial value and adopting a forward numerical integration method to obtain a corrected forward orbit initial value. And the reverse orbit initial value correction module is used for correcting the reverse orbit initial value of the orbit by using the tail epoch as an integral initial value and adopting a reverse numerical integration method to obtain a corrected reverse orbit initial value. And the forward reference orbit determination module is used for obtaining the forward reference orbit by taking the corrected initial value of the forward orbit as an integral initial value and adopting a forward numerical integration method. And the forward reference inter-satellite distance variability determining module is used for determining forward reference inter-satellite distance variability according to the forward reference orbit. And the reverse reference orbit determination module is used for obtaining a reverse reference orbit by taking the corrected reverse orbit initial value as an integral initial value and adopting a reverse numerical integration method. The reverse reference inter-satellite distance variability determining module is used for determining reverse reference inter-satellite distance variability according to the reverse reference orbit. The satellite gravity field inversion model building module is used for building a satellite gravity field inversion model according to the forward reference orbit, the forward reference inter-satellite distance variability, the backward reference orbit and the backward reference inter-satellite distance variability.

The satellite gravity field inversion model building module specifically comprises: the system comprises a forward observation equation building unit, a forward satellite gravity field inversion model determining unit, a reverse observation equation building unit, a reverse satellite gravity field inversion model determining unit and a satellite gravity field inversion model determining unit.

The forward observation equation building unit is used for building a forward observation equation according to the forward reference orbit and the forward reference inter-satellite distance variability. And the forward satellite gravity field inversion model determining unit is used for determining a forward satellite gravity field inversion model according to the forward observation equation. The reverse observation equation building unit is used for building a reverse observation equation according to the reverse reference orbit and the reverse reference inter-satellite distance variability. And the reverse satellite gravity field inversion model determining unit is used for determining a reverse satellite gravity field inversion model according to the reverse observation equation. The satellite gravity field inversion model determining unit is used for determining a satellite gravity field inversion model according to the forward satellite gravity field inversion model and the reverse satellite gravity field inversion model.

In order to acquire the accuracy and completeness of the obtained data, the system further comprises: the device comprises a data acquisition module and a preprocessing module.

The data acquisition module is used for acquiring the gravity satellite actual measurement data in the gravity field of the satellite to be inverted; the gravity satellite measured data comprises: orbit data, accelerometer data, star sensor data, and inter-star range data. And the preprocessing module is used for preprocessing the acquired actually measured data of the gravity satellite.

The embodiments in the present description are described in a progressive manner, each embodiment focuses on differences from other embodiments, and the same and similar parts among the embodiments are referred to each other. For the system disclosed by the embodiment, the description is relatively simple because the system corresponds to the method disclosed by the embodiment, and the relevant points can be referred to the method part for description.

The principles and embodiments of the present invention have been described herein using specific examples, which are provided only to help understand the method and the core concept of the present invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, the specific embodiments and the application range may be changed. In view of the above, the present disclosure should not be construed as limiting the invention.