CN111708095A - Satellite gravity field inversion method and system based on bidirectional integral - Google Patents

Satellite gravity field inversion method and system based on bidirectional integral Download PDF

Info

Publication number
CN111708095A
CN111708095A CN202010446920.1A CN202010446920A CN111708095A CN 111708095 A CN111708095 A CN 111708095A CN 202010446920 A CN202010446920 A CN 202010446920A CN 111708095 A CN111708095 A CN 111708095A
Authority
CN
China
Prior art keywords
satellite
orbit
gravity field
reverse
epoch
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010446920.1A
Other languages
Chinese (zh)
Other versions
CN111708095B (en
Inventor
肖云
王丽兵
成伟
赵玉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
61540 Troops of PLA
Original Assignee
61540 Troops of PLA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 61540 Troops of PLA filed Critical 61540 Troops of PLA
Priority to CN202010446920.1A priority Critical patent/CN111708095B/en
Publication of CN111708095A publication Critical patent/CN111708095A/en
Application granted granted Critical
Publication of CN111708095B publication Critical patent/CN111708095B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G01V7/02Details
    • G01V7/06Analysis or interpretation of gravimetric records

Abstract

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.

Description

Satellite gravity field inversion method and system based on bidirectional integral
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 as
Figure BDA0002506226330000081
Velocity vector of
Figure BDA0002506226330000082
Figure BDA0002506226330000083
Is a vector to be estimated, an
Figure BDA0002506226330000084
Wherein the content of the first and second substances,
Figure BDA0002506226330000085
x, y and z respectively represent three axes of a coordinate system and are three components of the position. v. ofx,vy,vzIs a velocity three component. bx,by,bzThe deviation is three components. k is a radical ofx,ky,kzIs a scale three component. Cnm、SnmAll gravity field model bit coefficients, n is the model order, m is the model number, when m is 0, coefficient Cn0Called the band harmonic coefficient, when n ═ m, coefficient Cnm、SnmCalled the fan harmonic coefficient, when n ≠ m, the coefficient Cnm、SnmReferred to as the field harmonic coefficients.
Then the corresponding state vector
Figure BDA0002506226330000091
Initial value of vector to be estimated
Figure BDA0002506226330000092
Partial derivative phi of, and
Figure BDA0002506226330000093
initial value of vector to be estimated
Figure BDA0002506226330000094
Partial derivatives of
Figure BDA0002506226330000095
Comprises the following steps:
Figure BDA0002506226330000096
Figure BDA0002506226330000097
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 values
Figure BDA0002506226330000098
Wherein
Figure BDA0002506226330000099
The integral vector is
Figure BDA00025062263300000910
And
Figure BDA00025062263300000911
respectively 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 t0The initial value of the integral vector is recorded as
Figure BDA00025062263300000912
And
Figure BDA00025062263300000913
then:
Figure BDA00025062263300000914
Figure BDA00025062263300000915
wherein the content of the first and second substances,
Figure BDA00025062263300000916
the right function is integrated for the start time,
Figure BDA00025062263300000917
is tkThe vector of the integral of the time of day,
Figure BDA00025062263300000918
is tkIntegral 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 t0The initial value of the integral vector is recorded as
Figure BDA00025062263300000919
And
Figure BDA00025062263300000920
then:
Figure BDA0002506226330000101
Figure BDA0002506226330000102
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:
Figure BDA0002506226330000103
Figure BDA0002506226330000104
Δriis the (i ═ A or B) orbital residual of the satellite i,
Figure BDA0002506226330000105
Is the inter-satellite distance-variability residual error,
Figure BDA0002506226330000106
is an inter-satellite distance-variability reference value, ri0,Bi0,Ki0Respectively, initial state parameters, accelerometer bias and scale parameters for satellite i (i ═ a or B), β are gravity field coefficients,
Figure BDA0002506226330000107
solved in the process of integrating the orbit motion equation,
Figure BDA0002506226330000108
the partial derivative of the inter-satellite distance variability to the initial state parameter and the gravity field coefficient is determined by
Figure BDA0002506226330000109
And (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:
Figure BDA0002506226330000111
wherein, ViFor observation errors, LiFor orbital residuals or inter-satellite range-rate residuals, AiThe 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:
Figure BDA0002506226330000112
wherein A isTIs 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:
Figure BDA0002506226330000113
xL、xGrespectively 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), N11,N22Are sub-matrices corresponding to local parameters and global parameters,
Figure BDA0002506226330000114
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
Figure BDA0002506226330000115
Similarly, each arc intersatellite speed measurement observation KBRR (K-band range-rate) can form a normal equation
Figure BDA0002506226330000121
In combination with the orbit data and the KBRR data, there is the following relationship:
Figure BDA0002506226330000122
wherein the content of the first and second substances,
Figure BDA0002506226330000123
Figure BDA0002506226330000124
σ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
Figure BDA0002506226330000125
The above formula is denoted as
Figure BDA0002506226330000126
Solving the equation of the reduced law to obtain the global parameter variable quantity of
Figure BDA0002506226330000127
Will be Δ xGBy substituting the following formula, the local parameter variation quantity Deltax that can be solvedL
Figure BDA0002506226330000128
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:
Figure BDA0002506226330000141
wherein the content of the first and second substances,
Figure BDA0002506226330000142
for post-fusion model bit coefficients, β1For the forward estimated bit coefficients, β2For the backward estimated bit coefficients, W1Estimating the weight of the bit coefficient for the forward direction, W2Is 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 ways1=W2=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 Cnm、SnmIn 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
Figure BDA0002506226330000151
Figure BDA0002506226330000161
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.

Claims (8)

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.
CN202010446920.1A 2020-05-25 2020-05-25 Satellite gravitational field inversion method and system based on bidirectional integration Active CN111708095B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010446920.1A CN111708095B (en) 2020-05-25 2020-05-25 Satellite gravitational field inversion method and system based on bidirectional integration

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010446920.1A CN111708095B (en) 2020-05-25 2020-05-25 Satellite gravitational field inversion method and system based on bidirectional integration

Publications (2)

Publication Number Publication Date
CN111708095A true CN111708095A (en) 2020-09-25
CN111708095B CN111708095B (en) 2023-07-28

Family

ID=72537340

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010446920.1A Active CN111708095B (en) 2020-05-25 2020-05-25 Satellite gravitational field inversion method and system based on bidirectional integration

Country Status (1)

Country Link
CN (1) CN111708095B (en)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102393535A (en) * 2011-07-20 2012-03-28 中国科学院测量与地球物理研究所 Satellite gravity inversion method based on two-star energy interpolation principle
CN103076640A (en) * 2013-01-17 2013-05-01 中国科学院测量与地球物理研究所 Method for inverting earth gravitational field by using variance-covariance diagonal tensor principle
WO2015042754A1 (en) * 2013-09-29 2015-04-02 清华大学 Method for parsing and calculating performance of satellite gravity field measurement by low-to-low satellite-to-satellite tracking
CN108267792A (en) * 2018-04-13 2018-07-10 武汉大学 Building global gravitational field model inversion method
CN108845366A (en) * 2018-06-04 2018-11-20 中国人民解放军61540部队 The system of harmonic analysis model and its modeling method of Satellite gravity field tensor diagonal line three-component inverting earth gravitational field

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102393535A (en) * 2011-07-20 2012-03-28 中国科学院测量与地球物理研究所 Satellite gravity inversion method based on two-star energy interpolation principle
CN103076640A (en) * 2013-01-17 2013-05-01 中国科学院测量与地球物理研究所 Method for inverting earth gravitational field by using variance-covariance diagonal tensor principle
WO2015042754A1 (en) * 2013-09-29 2015-04-02 清华大学 Method for parsing and calculating performance of satellite gravity field measurement by low-to-low satellite-to-satellite tracking
CN108267792A (en) * 2018-04-13 2018-07-10 武汉大学 Building global gravitational field model inversion method
CN108845366A (en) * 2018-06-04 2018-11-20 中国人民解放军61540部队 The system of harmonic analysis model and its modeling method of Satellite gravity field tensor diagonal line three-component inverting earth gravitational field

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
张兴福 等: "低轨卫星初始状态误差对反演地球重力场的影响" *

Also Published As

Publication number Publication date
CN111708095B (en) 2023-07-28

Similar Documents

Publication Publication Date Title
CN108225337B (en) Star sensor and gyroscope combined attitude determination method based on SR-UKF filtering
Pail et al. First GOCE gravity field models derived by three different approaches
CN107024212B (en) A kind of deep space probe X-ray pulsar/time difference astronomy doppler combined navigation method
Xu et al. Multiple parameter regularization: numerical solutions and applications to the determination of geopotential from precise satellite orbits
Pitjeva Modern numerical ephemerides of planets and the importance of ranging observations for their creation
CN104121927A (en) Inertial measurement unit calibration method applicable to low-accuracy no-azimuth-reference single-axis transposition equipment
CN113375677A (en) Spacecraft speed fixing method based on pulsar observation
CN112240941B (en) Relative calibration method and system for gravity satellite-borne accelerometer
CN103017787A (en) Initial alignment method suitable for rocking base
CN108959734A (en) One kind being based on real-time recursion solar light pressure torque discrimination method and system
CN105988129A (en) Scalar-estimation-algorithm-based INS/GNSS combined navigation method
CN103123487B (en) A kind of spacecraft attitude determination method
CN111580163A (en) Full waveform inversion method and system based on non-monotonic search technology
CN114595946A (en) Sea area theory minimum tide level calculation method, system, equipment and medium
CN103983274A (en) Inertial measurement unit calibration method suitable for low-precision no-azimuth reference biaxial transfer equipment
CN111708095A (en) Satellite gravity field inversion method and system based on bidirectional integral
CN107421533B (en) A kind of deep space probe X-ray pulsar TOA/DTOA Combinated navigation method
CN113449384A (en) Attitude determination method based on central error entropy criterion extended Kalman filtering
CN112632454A (en) MEMS gyro filtering method based on adaptive Kalman filtering algorithm
CN109557594B (en) Gravity reference graph time-varying correction method and system based on gravity abnormal time variation
CN103994775A (en) Inertia measurement unit calibration method for low precision dual-axis rotation-dwell equipment having orientation reference
CN108692727B (en) Strapdown inertial navigation system with nonlinear compensation filter
CN113405567B (en) Gravity satellite star sensor mounting matrix on-orbit calibration method and system
Leung et al. A comparison of the pseudo-linear and extended kalman filters for spacecraft attitude estimation
Gutsche et al. Podcast: Precise orbit determination software for leo satellites

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant