CN107065025B - A kind of orbital elements estimation method based on gravimetric field gradient invariant - Google Patents
A kind of orbital elements estimation method based on gravimetric field gradient invariant Download PDFInfo
- Publication number
- CN107065025B CN107065025B CN201710024671.5A CN201710024671A CN107065025B CN 107065025 B CN107065025 B CN 107065025B CN 201710024671 A CN201710024671 A CN 201710024671A CN 107065025 B CN107065025 B CN 107065025B
- Authority
- CN
- China
- Prior art keywords
- formula
- epoch
- orbit
- above formula
- follows
- 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.)
- Active
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V7/00—Measuring gravitational fields or waves; Gravimetric prospecting or detecting
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Navigation (AREA)
- Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
Abstract
A kind of orbital elements estimation method based on gravimetric field gradient invariant, its step are as follows: one: preparation；Two: decomposition and characteristic value of the ideal gravity gradient tensor under (East-North-Up, the ENU) coordinate system of northeast day；Three: solving decomposition and characteristic value of the measurement epoch J2 model gravity field tensor under ENU coordinate system；Four: using each measurement epoch r of J2 gravity gradient Searching Matrix Eigen Value and；Five: utilizing and calculate preliminary orbit element；Six: orbital elements are smooth.Gravity gradient matrix invariant information measured by satellite transit in the process has been used by what the above step present invention innovated, by iteratively solving out satellite to the geometric distance and reduced latitude in the earth's core, and give in this, as observed quantity using semi-major axis, eccentricity, orbit inclination angle, the perigee angle of depression and the true anomaly in orbital tracking as the measurement equation of state parameter；And the perturbation kinetics equation of orbital tracking is innovatively introduced, five orbital trackings of initial epoch are gone out using batch processing least-squares estimation；This method has the advantages that autonomy height, strong antijamming capability, construction cost are low etc., and the advantage that in deep-space detection field there is conventional method not have.
Description
Technical field
The present invention provides a kind of orbital elements estimation method based on gravimetric field gradient invariant, it is related to one kind to in-orbit
The method that satellite carries out orbital elements estimation using gravimetric field gradient measurement data, belongs to field of navigation technology.
Background technique
Earth gravitational field is the spatial dimension acted on by terrestrial gravitation, and irregular due to earth interior Mass Distribution
Property, so that earth gravitational field is not the field of force by a simple rule variation.In order to establish earth gravity field model, traditional object
Manage a series of gravitational field theory and methods of the observation data formation and development that geodesy is collected according to earth surface, such as weight
Power gesture, gravity bit function, geoid etc..The earth gravity field model formed due to the limitation of observation method, this stage
It is relatively rough, relatively accurate terrestrial gravitation field distribution can not be provided, so that the application using terrestrial gravitation field characteristic is limited.
Into the second half in 20th century, new measuring technique (such as electromagnetic distance measurement, artificial earth satellite positioning system and very long baseline are utilized
Interferometry (VLBI) etc.) and advanced gravity measurement instrument to establish more accurate earth gravity field model provide branch
It holds.Nowadays, due to the foundation of high-precision earth gravity field model, navigated using gravity gradient, determine appearance become it is numerous
The hot spot that country, scholar study.
Since the sixties in last century, there is scholar just to propose to utilize gravity gradiometer aided inertial navigation.Gravity gradient is auxiliary
Help inertial navigation, refer to using gravity gradient value correct inertial navigation system (Inertial Navigation System,
INS the position error) accumulated at any time has many advantages, such as independence strong, good concealment, is not limited by region and time domain.But by
Make gravitational field secondary navigation system not in the inherent shortcoming (such as gyroscopic drift and accelerometer bias are set) of inertia device
It can fundamentally solve the problems, such as the INS accumulation of error.Moreover, in the navigation research of these gravity gradients, main application scenarios limitation
In underwater submarine and aviation aircraft, gravity gradient navigation is without becoming a kind of independent air navigation aid.
In recent years, people start the gravity gradient navigation for being conceived to satellite in orbit, provide one kind for the autonomous operation of satellite
New departure.Traditional satellite navigation orbit determination method is mainly based upon the orbit determination method of earth station's telemetering.However, using earth station pair
Satellite carries out orbit determination, and there are many restrictions, such as construction cost, ground station failure, data-handling capacity, data movement capacity, region
Limitation, political factor etc..For the China for possessing a large amount of launch missions at present, the quantity and processing capacity of earth station are simultaneously
Current needs cannot be fully met.It is therefore desirable to propose the novel orbit determination method of one kind as the benefit of traditional orbit determination method
Fill, gravity gradient navigation becomes one of the effective scheme for avoiding the above drawback, it is low with construction cost, do not depend on completely it is outer
The information that boundary provides, independence is strong on star, not the advantage vulnerable to extraneous factor interference.
Currently, the spacefaring nations such as the U.S., Russia and China are proposed the survey plan for deep layer space.Work as spy
After device is surveyed far from the earth, continue to will appear using traditional earth station's telemetering orbit determination that signal strength is faint, and delay is serious, by various
The problems such as interference, it will bring huge difficulty to deep space probe orbit determination in real time.And if row will be detected in advance by having had been established
The gravity field model of star, then gravimetric field gradient navigation will provide autonomous, real-time, reliable navigation for satellite, orbit determination is tied
Fruit.
To sum up, due to the foundation of current high-precision earth gravity field model and the appearance of advanced gravity measurement instrument
Support is provided for gravity gradient navigation.With going deep into for research, is become using gravity gradient for spacecraft progress orbit determination and worked as
The hot spot of preceding research.And since it allows its research to be provided with necessity compared to advantage possessed by traditional orbit determination method and meaning
Property.The present invention just proposes a kind of orbital elements estimation method based on gravimetric field gradient invariant, uses gravimetric field gradient square
Invariant information in battle array estimates five elements in six element of spacecraft orbit.
Summary of the invention
(1) goal of the invention:
The present invention is intended to provide a kind of novel autonomous orbit determination, air navigation aid, this method innovation has used satellite transit
Measured gravity gradient matrix invariant information in the process, by iteratively solving out geometric distance and the earth's core of the satellite to the earth's core
Latitude, and in this, as observed quantity give in orbital tracking semi-major axis, eccentricity, orbit inclination angle, the perigee angle of depression and
True anomaly is the measurement equation of state parameter.In order to improve the estimated accuracy of orbital tracking, innovation introduces orbital tracking
Perturbation kinetics equation, five orbital trackings of initial epoch are gone out using batch processing least-squares estimation.This method has certainly
The advantages such as main degree height, strong antijamming capability, construction cost be low, and in deep-space detection field there is conventional method not had
Advantage.The invention proposes a complete orbital tracking estimation method processes, and for the difficulty in computation during reducing,
Give relatively simple calculation method.
(2) technical solution
A kind of orbital elements estimation method based on gravimetric field gradient invariant of the present invention, its step are as follows.
Step 1: preparation
Terrestrial gravitation potential function
In the Primary Study spacecraft characteristics of motion, the earth is regarded as ideal sphere, the gravitation that it generates spacecraft
It is directed toward the earth's core, size F is directly proportional to the quality m of spacecraft, and square is inversely proportional with the earth's core to spacecraft distance r:
In above formula, F is terrestrial gravitation suffered by spacecraft, and G is terrestrial gravitation constant, and m is spacecraft mass, r be the earth's core extremely
The geometric distance of spacecraft.F is exactly gravitational acceleration g divided by m.Two constants are merged μ=GM, then obtain following formula:
In above formula, μ becomes Gravitational coefficient of the Earth.
Earth gravitational field becomes center gravitational field in the case of this.Gravitational acceleration is expressed as the gradient of gravitational potential, then in
The gravitational potential function of heart gravitational field is:
The actual earth be not it is spherically symmetric, there is degree of flattening, pyriform and equator deformation etc., thus gravitational acceleration is
Distance r, latitude (reduced latitude)With the function of longitude λ, and r,λ is connected coordinate system S in the eartheThe spherical coordinates of middle position.
In order to strictly and easily describe the distribution of normal gravity, it is expressed as the gradient of terrestrial gravitation potential function U:
Grad (*) is to seek gradient operator in above formula.
It is given below, is only considering the aspherical J of the earth2Gravitational field expression formula under the influence of:
R=6378.1363km in above formula is terrestrial equator radius.
Kepler orbit elements
Spacecraft orbit is observed in inertial space, needs to define orbital elements.Orbital elements, also known as orbital tracking are to determine
Determine one group of constant of the motion feature of Keplerian orbit.After considering orbit perturbation, orbital tracking is no longer constant, be can be used as
State variable.
Orbital tracking one shares 6, indicates and meaning is as follows:
A: elliptic orbit half-court axis；
E: elliptic orbit eccentricity；
Ω: right ascension of ascending node；The point that spacecraft orbit passes through equator from south to north is known as ascending node B.In equatorial plane, by
The angle that the first point of Aries goes to eastwards ascending node B is known as right ascension of ascending node, and effective range is 0 to 360 degree；
I: orbit inclination angle；The angle of orbit plane and equatorial plane.Effective range is 0 to 180 degree；
tp: near-earth point moment；Spacecraft can also use moment in epoch t by the perigean moment0True anomaly θ
(t0) or eccentric anomaly E (t0) or mean anomaly M (t0)；
ω: argument of perigee；In orbit plane, there is ascending node to go to perigean angle.Effective range is 0 to 360 degree.
Their effect is respectively: semi-major axis of orbit a and eccentric ratio e determine the size and shape of track；Right ascension of ascending node
Ω and orbit inclination angle i determines orbit plane in the orientation of inertial space；Argument of perigee ω determines the line of apsides in orbit plane
Position；Mean anomaly M (the t at moment in epoch0) determine the starting point of time.
Orbit perturbation equation
The track of spacecraft is not Keplerian orbit under perturbative force effect, and orbital elements change over time: a
(t), e (t), Ω (t), i (t) ....
Orbit perturbation equation is given below:
In above formula, there are identical meanings with the symbol occurred above；frIndicate perturbation acceleration in geocentric orbit coordinate
The radial component of system；fuIndicate perturbation acceleration in the cross stream component of geocentric orbital reference system；fhIndicate perturbation acceleration in the earth's core
The secondary normal component of orbital coordinate system；P is semi-latus rectum.
Step 2: decomposition and spy of the ideal gravity gradient tensor under (East-North-Up, the ENU) coordinate system of northeast day
Value indicative
Shown in ideal sphere gravity field model such as formula (3), r has following expression in formula:
In above formula, x, y, z are three component projections that satellite position vector r is connected on coordinate system in the earth respectively.
Gravity gradient tensor can regard the gradient of gravitational acceleration as, and gravitational acceleration is the gradient of gravitational field again.Cause
And it is directly decomposed under ENU coordinate system:
In above formula, subscript E, N, U indicate the different elements that Γ is decomposed under ENU coordinate system.
Three characteristic values for obtaining gravity gradient matrix that can be convenient by formula (8):
It can know from formula (9), three characteristic values are linearly related two-by-two, only one characteristic quantity.
Step 3: decomposition and characteristic value of the measurement epoch J2 model gravity field tensor under ENU coordinate system are solved
J2 model gravity field tensor is provided by formula (5).It is sat using asking gradient operator to solve J2 model gravity field tensor in ENU
Decomposition under mark system:
Wherein, each component form is as follows:
By formula (10) solution matrixEigenvalue λ1, λ2, λ3, characteristic value withRelationship it is as follows:
λ can be verified1, λ2, λ3There is following relationship:
λ1+λ2+λ3=0 (13)
Therefore, only there are two independent quantities in three characteristic values.
Step 4: using each measurement epoch r of J2 gravity gradient Searching Matrix Eigen Value and
1) r initial value is calculated
Use the gravity gradient tensor matrix exgenvalue (12) of J2 model solution, and ideal sphere gravity field model characteristic value
Expression formula (9) calculates r initial value.
Eigenvalue λ in formula (12)1And λ2Optional one all can, bring formula (9) into, the anti-r that solves then has following initial value expression formula:
Or
2) it calculatesInitial value
Following calculating is done by formula (12):
It is available from formula (16)Expression formula:
It is available from formula (17)Another expression formula:
More thanTwo kinds of expression formulas, optional one.Exist for the sign of cosine value in formula (18) and formula (19)
This is not discussed.
3) it utilizesUpdate r
Utilize the λ in formula (12)2Expression formula updates r, by r0WithBring anti-solution r into3, expression formula is as follows:
4) r is utilized1It updates
Method one:
Following calculating is done with three characteristic values in formula (12):
By r1WithBring above formula (23) into, it is anti-to solve in above formula, it is as follows:
Method two:
Following calculating is done using formula (21) and formula (22):
By r1WithBring above formula (25) into, it is anti-to solve in above formula, it is as follows:
Step 5: using r andCalculate preliminary orbit element
1) semi-major axis and eccentricity are calculated
Maximum value r is searched in all result r for calculating epochmaxWith minimum value rmin, then there is following relationship:
ra=rmax, rp=rmin (27)
In above formula, subscript a indicates apogee, and subscript p indicates perigee.
Then semi-major axis:
Eccentricity:
2) orbit inclination angle is calculated
It is assumed herein that obtained by step 4The sign of value is learnt by other methods, in all meters
Calculate the result of epochMiddle search maximum value,
Orbit inclination angle:
3) argument of perigee and initial true anomaly are calculated
The r searched in 2)minIts corresponding latitude is denoted asThen argument of perigee ω is calculated by following formula:
Initial epoch true anomaly θ0, it is calculated by following formula:
In above formula,For first epoch calculated latitude.
Step 6: orbital elements are smooth
It is as follows using five orbital elements of initial epoch as parameter to be estimated:
X=[a0 ex0 ey0 i0 u0] (33)
In above formula, subscript " 0 " indicates initial epoch parameter.ex, ey, u expression formula is as follows:
The initial value of parameter x to be estimated, is provided by step 5.
Dynamics of orbits model is described by formula (6), semi-latus rectumThe perturbative force due to caused by J2 adds
Decomposition of the speed under geocentric orbital reference system is as follows:
Convolution (6) and formula (35), and the initial value provided by step 5 just constitute the differential equation group of known initial value
Solve problems can then acquire five orbital elements of any k epoch-making moment, x by numerical integration methodk=[ak exk eyk ik
uk], then the observational equation of k epoch-making moment is as follows:
In above formula,For the geometric distance and latitude of each epoch satellite obtained in step 4 to the earth's core,Table
Show that calculating error, subscript " k " indicate k-th of epoch.
Simply assumeIt is a Gaussian noise, then k epoch, it is as follows to measure covariance matrix:
In above formula, σλFor gravity gradient invariant error, error size will be according to the measurement noise and J2 of measuring instrument
Model error size is chosen.
The measured value z of each epoch is put together into the form of being written as follow:
In above formula, N is the number of measurement epoch.Each epoch measurement covariance matrix is write as following form:
Estimated state parameter is filtered using quasi- Sigma's least square batch processing, estimator is as follows:
In above formula, subscript j indicates the number of iterations,It indicates to utilizeThe calculated value of calculated measurement amount z with
Sensitive matrix.
In formula (40), H gusts of difficulty of solution are larger, and Z is write as about x=[a0 ex0 ry0 i0 u0] implicit function expression,
It is as follows:
Z=f (a0 ex0 ey0 i0 u0) (41)
Local derviation is asked to calculate sensitive matrix H by general analytical method, calculation amount is very huge.It therefore can be by the following method
It solves.
It is finished in iteration j, obtains state estimation valueThen carry out such as down-sampling:
In above formula, L is the dimension of x, P0For priori covariance matrix,Indicate P0The i-th column square root.Priori association
Variance matrix is given by:
Formula (42) are integrated by formula (6) and formula (36) transmitting, calculating measured value are as follows:
χI, j→γI, j, i=0 ..., L (44)
Then density matrix H can be obtained approximately through following formula:
In above formula,Indicate P0The square root of diagonal entry.
Through the above steps, a kind of orbital elements estimation method based on gravimetric field gradient invariant is proposed.This method
Using gravity gradient information measured by gravity gradiometer, and utilize ideal sphere gravity field model and the gravity gradient with J2
The characteristic value of model solution gravity gradient tensor gives the method for solving five orbital elements using this feature value, in order to mention
High solving precision introduces orbit perturbation kinetics equation and carrys out smooth five orbital elements.For computationally intensive part, give
The corresponding method for reducing difficulty in computation.To sum up, this method is completely independent of the input of external information, has very strong autonomous
Property, construction cost is low to have very strong application prospect in independent navigation and deep hole field of detecting.If the higher precision introduced
Gravitational field will further increase orbit determination accuracy, have very strong potentiality value.
(3) advantage
A kind of the advantages of orbital elements estimation method method based on gravimetric field gradient invariant provided by the invention, is:
1. innovation has used invariant information in gradiometry matrix, gives in conjunction with gravity field model in the present invention
The method for calculating orbital elements using invariant information is gone out.Determine that orbital elements do not need to additionally introduce other external worlds in the process
Information has very strong independence, anti-interference ability；
2. method proposed by the present invention, with the low advantage of cost of implementation.A large amount of money is needed different from earth station's orbit determination
Gold investment carries out system Construction, and this method need to only carry a gradiometry instrument on satellite can realize autonomous orbit determination,
Cost of implementation is low；
3. method proposed by the present invention has huge application prospect in deep-space detection field.Grasping planet to be surveyed
After accurate gravity field model, real-time independent navigation, orbit determination can be realized.
Detailed description of the invention
Fig. 1 is the method for the invention flow chart.
Fig. 2 be in the present invention using each measurement epoch r of J2 gravity gradient Searching Matrix Eigen Value andFlow chart
Specific embodiment
Specific implementation process of the invention is done further specifically below in conjunction with attached drawing 1, Fig. 2 and technical solution
It is bright.
Step 1: decomposition and characteristic value of the ideal gravity gradient tensor under (ENU) coordinate system of northeast day
Shown in ideal sphere gravity field model such as formula (3), r has following expression in formula:
In above formula, x, y, z are three component projections that satellite position vector r is connected on coordinate system in the earth respectively.
Gravity gradient tensor can regard the gradient of gravitational acceleration as, and gravitational acceleration is the gradient of gravitational field again.Cause
And it is directly decomposed under ENU coordinate system:
In above formula, subscript E, N, U indicate the different elements that Γ is decomposed under ENU coordinate system.
Three characteristic values for obtaining gravity gradient matrix that can be convenient by formula (47):
It can know from formula (48), three characteristic values are linearly related two-by-two, only one characteristic quantity.
Second box in the step respective figure 1.
Step 2: decomposition and characteristic value of the measurement epoch J2 model gravity field tensor under ENU coordinate system are solved
J2 model gravity field tensor is provided by formula (5).J2 model gravity field tensor is solved in ENU seat by asking gradient to calculate
Decomposition under mark system, obtains following formula:
Wherein, each component form is as follows:
By formula (49) solution matrix eigenvalue λ1, λ2, λ3, characteristic value withRelationship it is as follows:
λ can be verified1, λ2, λ3There is following relationship:
λ1+λ2+λ3=0 (52)
Therefore, only there are two independent quantities in three characteristic values.
Second box in the step respective figure 1.
Step 3: using each measurement epoch r of J2 gravity gradient Searching Matrix Eigen Value and
1) r initial value is calculated
Use the gravity gradient tensor matrix exgenvalue (51) of J2 model solution, and ideal sphere gravity field model characteristic value
Expression formula (48) calculates r initial value.
Eigenvalue λ in formula (51)1And λ2Optional one all can, bring formula (48) into, the anti-r that solves then has following initial value expression formula:
Or
Second box in the step respective figure 2.
2) it calculatesInitial value
Following calculating is done by formula (51):
It is available from formula (55)Expression formula:
Or, available from formula (56)Another expression formula:
More thanTwo kinds of expression formulas, optional one.Exist for the sign of cosine value in formula (57) and formula (58)
This is not discussed.
Third box in the step respective figure 2.
3) it utilizesUpdate r
Utilize the λ in formula (51)2Expression formula updates r, by λ0WithBring anti-solution r into3, expression formula is as follows:
4th box in the step respective figure 2.
4) r is utilized1It updates
Method one:
Following calculating is done with three characteristic values in formula (51):
By r1WithBring above formula (62) into, it is anti-to solve in above formula, it is as follows:
Method two:
Following calculating is done using formula (60) and formula (61):
By r1WithBring above formula (64) into, it is anti-to solve in above formula, it is as follows:
5th box in the step respective figure 2.
5) judge geometric distance r and latitudeWhether precision prescribed is reached
Circulation is jumped out if reaching and meeting precision, is executed step 4, is otherwise executed 3).
The 6th box in the step respective figure 2.
The 4th box in step 3 respective figure 1.
Step 4: judge whether measurement epoch calculates and finish
If all epoch, which calculate to finish, thens follow the steps five, step 2 is otherwise executed, calculates next epoch gravity gradient
Moment matrix characteristic value.
The 5th box in the step respective figure 1.
Step 5: using r andCalculate preliminary orbit element
1) semi-major axis and eccentricity are calculated
Maximum value r is searched in all result r for calculating epochmaxWith minimum value rmin, then there is following relationship:
ra=rmax, rp=rmin (66)
In above formula, subscript a indicates apogee, and subscript p indicates perigee.
Then semi-major axis:
Eccentricity:
2) orbit inclination angle is calculated
It is assumed herein that obtained by step 4The sign of value is learnt by other methods, in all meters
Calculate the result of epochMiddle search maximum value,
Orbit inclination angle:
3) argument of perigee and initial true anomaly are calculated
The r searched in 2)minIts corresponding latitude is denoted asThen argument of perigee ω is calculated by following formula:
Initial epoch true anomaly θ0, it is calculated by following formula:
In above formula,For first epoch calculated latitude.
Step 6: orbital elements are smooth
It is as follows using five orbital elements of initial epoch as parameter to be estimated:
X=[a0 ex0 ey0 i0 u0] (72)
In above formula, subscript " 0 " indicates initial epoch parameter.ex, ey, u expression formula is as follows:
The initial value of parameter x to be estimated, is provided by step 5.
Dynamics of orbits model is described by formula (6), semi-latus rectumThe perturbative force due to caused by J2 adds
Decomposition of the speed under geocentric orbital reference system is as follows:
Convolution (6) and formula (74), and the initial value provided by step 5 just constitute the differential equation group of known initial value
Solve problems can then acquire five orbital elements of any k epoch-making moment, x by numerical integration methodk=[ak exk eyk ik
uk], then the observational equation of k epoch-making moment is as follows:
In above formula,For the geometric distance and latitude of each epoch satellite obtained in step 3 to the earth's core,Table
Show that calculating error, subscript " k " indicate k-th of epoch.
Simply assumeIt is a Gaussian noise, then k epoch, it is as follows to measure covariance matrix:
In above formula, σλFor gravity gradient invariant error, error size will be according to the measurement noise and J2 of measuring instrument
Model error size is chosen.
The measured value z of each epoch is put together into the form of being written as follow:
In above formula, N is the number of measurement epoch.Each epoch measurement covariance matrix is write as following form:
Estimated state parameter is filtered using quasi- Sigma's least square batch processing, estimator is as follows:
In above formula, subscript j indicates the number of iterations,It indicates to utilizeThe calculated value of calculated measurement amount z with
Sensitive matrix.
In formula (79), H gusts of difficulty of solution are larger, and Z is write as about x=[a0 ex0 ey0 i0 u0] implicit function expression,
It is as follows:
Z=f (a0 ex0 ey0 i0 u0) (80)
Local derviation is asked to calculate sensitive matrix H by general analytical method, calculation amount is very huge.It therefore can be by the following method
It solves.
It is finished in iteration j, obtains state estimation valueThen carry out such as down-sampling:
In above formula, L is the dimension of x, P0For priori covariance matrix,Indicate P0The i-th column square root.Priori association side
Poor matrix is given by:
Formula (81) are integrated by formula (6) and formula (75) transmitting, calculating measured value are as follows:
χI, j→γI, j, i=0 ..., L (83)
Then density matrix H can be obtained approximately through following formula:
In above formula,Indicate P0The square root of diagonal entry.
Terminate the final smooth orbital elements result of iteration output when least-squares iteration reaches required precision.
Through the above steps, a kind of complete stream of orbital elements estimation method based on gravimetric field gradient invariant is given
Journey.The process uses gravity gradient matrix invariant information, introduces gravity field model and orbit perturbation equation, obtains five tracks
Element, and solution difficulty in computation is simplified, reduce the performance requirement to computer.To sum up, this method process independent of
The input of external information, in the requirement for meeting the spacecrafts autonomous operations such as current and later satellite, construction cost is far below ground
The construction of face station, and there is bright application prospect in following deep space exploration, independent navigation field.
Claims (1)
1. a kind of orbital elements estimation method based on gravimetric field gradient invariant, it is characterised in that: its step are as follows:
Step 1: preparation
1) terrestrial gravitation potential function
In the Primary Study spacecraft characteristics of motion, the earth is regarded as ideal sphere, it is directed toward the gravitation that spacecraft generates
The earth's core, size F is directly proportional to the quality m of spacecraft, and square is inversely proportional with the earth's core to spacecraft distance r:
In above formula, F is terrestrial gravitation suffered by spacecraft, and G is terrestrial gravitation constant, and m is spacecraft mass, and r is the earth's core to space flight
The geometric distance of device, F are exactly gravitational acceleration g divided by m, and two constants are merged μ=GM, then obtain following formula:
In above formula, μ becomes Gravitational coefficient of the Earth；
Earth gravitational field becomes center gravitational field in the case of this, gravitational acceleration is expressed as the gradient of gravitational potential, then center is drawn
The gravitational potential function in the field of force is:
The actual earth be not it is spherically symmetric, have flatness, pyriform and equator deformation, thus gravitational acceleration be distance r,
LatitudeWith the function of longitude λ, and r,λ is connected coordinate system S in the eartheThe spherical coordinates of middle position, in order to strictly and
The distribution for easily describing normal gravity is expressed as it the gradient of terrestrial gravitation potential function U:
Grad (*) is to seek gradient operator in above formula；
It is given below, is only considering the aspherical J of the earth2Gravitational field expression formula under the influence of:
Re=6378.1363km in above formula is terrestrial equator radius；
2) Kepler orbit elements
Spacecraft orbit is observed in inertial space, needs to define orbital elements；Orbital elements, also known as orbital tracking are to determine to open
One group of constant of the general motion feature for strangling track；After considering orbit perturbation, orbital tracking is no longer constant, can be used as state
Variable；
Orbital tracking one shares 6, indicates and meaning is as follows:
A: elliptic orbit half-court axis；
E: elliptic orbit eccentricity；
Ω: right ascension of ascending node；The point that spacecraft orbit passes through equator from south to north is known as ascending node B；In equatorial plane, by the Spring Equinox
The angle that point goes to eastwards ascending node B is known as right ascension of ascending node, and effective range is 0 to 360 degree；
I: orbit inclination angle；The angle of orbit plane and equatorial plane, effective range are 0 to 180 degree；
tp: near-earth point moment；Spacecraft passes through the perigean moment；
ω: argument of perigee；In orbit plane, there is ascending node to go to perigean angle, effective range is 0 to 360 degree；
Their effect is respectively: semi-major axis of orbit a and eccentric ratio e determine the size and shape of track；Right ascension of ascending node Ω and
Orbit inclination angle i determines orbit plane in the orientation of inertial space；Argument of perigee ω determines position of the line of apsides in orbit plane；
Mean anomaly M (the t at moment in epoch0) determine the starting point of time；
3) orbit perturbation equation
The track of spacecraft is not Keplerian orbit under perturbative force effect, and orbital elements change over time: a (t), e
(t), Ω (t), i (t), tp, ω (t)；
Orbit perturbation equation is given below:
In above formula, there are identical meanings with the symbol occurred above；frIndicate perturbation acceleration in the diameter of geocentric orbital reference system
To component；fuIndicate perturbation acceleration in the cross stream component of geocentric orbital reference system；fhIndicate that perturbation acceleration is sat in geocentric orbit
Mark the secondary normal component of system；P is semi-latus rectum；
Step 2: decomposition and characteristic value of the ideal gravity gradient tensor under (East-North-Up, the ENU) coordinate system of northeast day
Shown in ideal sphere gravity field model such as formula (3), r has following expression in formula:
In above formula, x, y, z are three component projections that satellite position vector r is connected on coordinate system in the earth respectively；
Gravity gradient tensor can regard the gradient of gravitational acceleration as, and gravitational acceleration is the gradient of gravitational field again, thus straight
It connects and is decomposed under ENU coordinate system:
In above formula, subscript E, N, U are indicatedThe different elements decomposed under ENU coordinate system；
Three characteristic values for obtaining gravity gradient matrix that can be convenient by formula (8):
It can know from formula (9), three characteristic values are linearly related two-by-two, only one characteristic quantity；
Step 3: decomposition and characteristic value of the measurement epoch J2 model gravity field tensor under ENU coordinate system are solved
J2 model gravity field tensor is provided by formula (5), using ask gradient operator solve J2 model gravity field tensor in ENU coordinate system
Under decomposition:
Wherein, each component form is as follows:
By formula (10) solution matrixEigenvalue λ1,λ2,λ3, characteristic value and r,Relationship it is as follows:
λ can be verified1, λ2, λ3There is following relationship:
λ1+λ2+λ3=0 (13)
Therefore, only there are two independent quantities in three characteristic values；
Step 4: using each measurement epoch r of J2 gravity gradient Searching Matrix Eigen Value and
1) r initial value is calculated
It is reached using the gravity gradient tensor matrix exgenvalue (12) of J2 model solution, and the ideal sphere gravity field model list of feature values
Formula (9) calculates r initial value；
Eigenvalue λ in formula (12)1And λ2Optional one all can, bring formula (9) into, the anti-r that solves then has following initial value expression formula:
Or
2) it calculatesInitial value
Following calculating is done by formula (12):
It is available from formula (16)Expression formula:
It is available from formula (17)Another expression formula:
More thanTwo kinds of expression formulas, optional one；
3) it utilizesUpdate r
Utilize the λ in formula (12)2Expression formula updates r, by r0WithBring anti-solution r into3, expression formula is as follows:
4) r is utilized1It updates
Method one:
Following calculating is done with three characteristic values in formula (12):
By r1WithBring above formula (23) into, it is anti-to solve in above formula, it is as follows:
Method two:
Following calculating is done using formula (21) and formula (22):
By r1WithBring above formula (25) into, it is anti-to solve in above formula, it is as follows:
Step 5: using r andCalculate preliminary orbit element
1) semi-major axis and eccentricity are calculated
Maximum value r is searched in all result r for calculating epochmaxWith minimum value rmin, then there is following relationship:
ra=rmax,rp=rmin (27)
In above formula, subscript a indicates apogee, and subscript p indicates perigee；
Then semi-major axis:
Eccentricity:
2) orbit inclination angle is calculated
It is assumed herein that obtained by step 4The sign of value has been learnt, in all results for calculating epochIn
Maximum value is searched for,
Orbit inclination angle:
3) argument of perigee and initial true anomaly are calculated
The r searched in 2)minIts corresponding latitude is denoted asThen argument of perigee ω is calculated by following formula:
Initial epoch true anomaly θ0, it is calculated by following formula:
In above formula,For first epoch calculated latitude；
Step 6: orbital elements are smooth
It is as follows using five orbital elements of initial epoch as parameter to be estimated:
X=[a0 ex0 ey0 i0 u0] (33)
In above formula, subscript " 0 " indicates initial epoch parameter, ex,ey, u expression formula is as follows:
ex=ecos ω
ey=esin ω (34)
U=ω+θ
The initial value of parameter x to be estimated, is provided by step 5；
Dynamics of orbits model is described by formula (6), semi-latus rectumThe perturbative force acceleration due to caused by J2
Decomposition under geocentric orbital reference system is as follows:
Convolution (6) and formula (35), and the initial value provided by step 5 just constitute the Systems of Ordinary Differential Equations of known initial value
Problem can then acquire five orbital elements of any k epoch-making moment, x by numerical integration methodk=[ak exk eyk ik uk],
Then the observational equation of k epoch-making moment is as follows:
In above formula,Geometric distance and latitude for each epoch satellite obtained in step 4 to the earth's core, Δ r,It indicates
Error is calculated, subscript " k " indicates k-th of epoch；
Simply assume Δ r,It is a Gaussian noise, then k epoch, it is as follows to measure covariance matrix:
In above formula, σλFor gravity gradient invariant error, error size will be missed according to the measurement noise and J2 model of measuring instrument
Poor size is chosen；
The measured value z of each epoch is put together into the form of being written as follow:
In above formula, N is the number of measurement epoch；Each epoch measurement covariance matrix is write as following form:
Estimated state parameter is filtered using quasi- Sigma's least square batch processing, estimator is as follows:
In above formula, subscript j indicates the number of iterations,HjIt indicates to utilizeThe calculated value and sensitivity square of calculated measurement amount z
Battle array；
In formula (40), H gusts of difficulty of solution are larger, and Z is write as about x=[a0 ex0 ey0 i0 u0] implicit function expression, it is as follows:
Z=f (a0 ex0 ey0 i0 u0) (41)
Local derviation is asked to calculate sensitive matrix H by general analytical method, calculation amount is very huge, therefore can ask by the following method
Solution；
It is finished in iteration j, obtains state estimation valueThen carry out such as down-sampling:
In above formula, L is the dimension of x, P0For priori covariance matrix,Indicate P0The i-th column square root；Priori covariance square
Battle array is given by:
Formula (42) are integrated by formula (6) and formula (36) transmitting, calculating measured value are as follows:
χi,j→γi,j, i=0 ..., L (44)
Then sensitive matrix H can be obtained approximately through following formula:
In above formula,Indicate P0The square root of i-th of diagonal entry, i=1,2,3,4,5.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710024671.5A CN107065025B (en) | 2017-01-13 | 2017-01-13 | A kind of orbital elements estimation method based on gravimetric field gradient invariant |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710024671.5A CN107065025B (en) | 2017-01-13 | 2017-01-13 | A kind of orbital elements estimation method based on gravimetric field gradient invariant |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107065025A CN107065025A (en) | 2017-08-18 |
| CN107065025B true CN107065025B (en) | 2019-04-23 |
Family
ID=59599280
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710024671.5A Active CN107065025B (en) | 2017-01-13 | 2017-01-13 | A kind of orbital elements estimation method based on gravimetric field gradient invariant |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107065025B (en) |
Families Citing this family (17)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108415244B (en) * | 2017-12-28 | 2020-12-18 | 中国空间技术研究院 | Multi-degree-of-freedom electrostatic suspension system geometric symmetry approximation method |
| CN108548542B (en) * | 2018-07-13 | 2021-09-28 | 北京航空航天大学 | Near-earth orbit determination method based on atmospheric resistance acceleration measurement |
| CN108919370B (en) * | 2018-07-25 | 2019-11-29 | 赛德雷特(珠海)航天科技有限公司 | A kind of positioning device and method based on gravitation field measurement |
| CN109738919B (en) * | 2019-02-28 | 2020-12-15 | 西安开阳微电子有限公司 | Method for autonomous ephemeris prediction of GPS receiver |
| CN110053788B (en) * | 2019-03-15 | 2022-05-13 | 中国西安卫星测控中心 | Constellation long-term retention control frequency estimation method considering complex perturbation |
| CN110378012B (en) * | 2019-07-16 | 2021-07-16 | 上海交通大学 | Strict regression orbit design method, system and medium considering high-order gravity field |
| CN110442831B (en) * | 2019-07-31 | 2023-03-24 | 中国人民解放军国防科技大学 | Space non-cooperative target space-based search method based on nonlinear deviation evolution |
| CN110705022B (en) * | 2019-08-30 | 2021-11-02 | 中国矿业大学 | Sparse spherical radial basis function local gravity field modeling method |
| CN110967041B (en) * | 2019-12-18 | 2021-09-14 | 自然资源部国土卫星遥感应用中心 | Tensor invariant theory-based satellite gravity gradient data precision verification method |
| CN111505679B (en) * | 2020-04-20 | 2022-05-03 | 中国科学院国家空间科学中心 | Satellite-borne GNSS-based LEO initial orbit determination method |
| CN111721303B (en) * | 2020-07-01 | 2022-09-13 | 中国人民解放军陆军装甲兵学院 | Spacecraft inertial navigation method, system, medium and equipment based on gravitational field |
| CN112762924B (en) * | 2020-12-23 | 2023-07-14 | 北京机电工程研究所 | Navigation positioning method based on gravity gradient-topography heterologous data matching |
| CN113433596B (en) * | 2021-06-25 | 2022-09-16 | 中国船舶重工集团公司第七0七研究所 | Gravity gradient dynamic measurement filtering method based on spatial domain |
| CN115327653B (en) * | 2022-08-15 | 2023-04-04 | 自然资源部国土卫星遥感应用中心 | Tensor invariant theory-based satellite gravity gradient gross error detection method |
| CN115792982B (en) * | 2022-11-07 | 2023-10-20 | 北京航空航天大学合肥创新研究院(北京航空航天大学合肥研究生院) | Beidou navigation satellite broadcast ephemeris parameter fitting method and storage medium |
| CN115808881B (en) * | 2023-01-21 | 2023-04-21 | 中国科学院数学与***科学研究院 | Method for estimating on-orbit quality of non-trailing satellite and self-adaptive control method |
| CN116108319B (en) * | 2023-04-10 | 2023-08-11 | 中国人民解放军32035部队 | Orbit forecasting method for constant thrust mode continuous maneuvering satellite |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104061932A (en) * | 2014-06-10 | 2014-09-24 | 中国空间技术研究院 | Method for navigation positioning by using gravitation vector and gradient tensor |
| CN105138005A (en) * | 2015-06-16 | 2015-12-09 | 北京航空航天大学 | Method for determining relative orbit elements based on inter-satellite distance |
| CN106092099A (en) * | 2016-06-02 | 2016-11-09 | 哈尔滨工业大学 | Spacecraft is relative to positional increment orbit determination method |
| CN106202640A (en) * | 2016-06-28 | 2016-12-07 | 西北工业大学 | Day ground three body gravitational fields in halo orbit spacecraft bias track method for designing |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8359920B2 (en) * | 2009-05-15 | 2013-01-29 | Lockheed Martin Corp. | Gravity sensing instrument |
-
2017
- 2017-01-13 CN CN201710024671.5A patent/CN107065025B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104061932A (en) * | 2014-06-10 | 2014-09-24 | 中国空间技术研究院 | Method for navigation positioning by using gravitation vector and gradient tensor |
| CN105138005A (en) * | 2015-06-16 | 2015-12-09 | 北京航空航天大学 | Method for determining relative orbit elements based on inter-satellite distance |
| CN106092099A (en) * | 2016-06-02 | 2016-11-09 | 哈尔滨工业大学 | Spacecraft is relative to positional increment orbit determination method |
| CN106202640A (en) * | 2016-06-28 | 2016-12-07 | 西北工业大学 | Day ground three body gravitational fields in halo orbit spacecraft bias track method for designing |
Non-Patent Citations (3)
| Title |
|---|
| Gravity Gradient Tensor Eigendecomposition for Spacecraft Positioning;Pei Chen 等;《Journal of guidance，control，and dynamics》;20151130;第38卷(第11期);第2200-2206页 |
| 利用引力梯度不变量解算的GOCE引力场模型;于锦海 等;《中国科学：地球科学》;20121231;第42卷(第9期);第1450-1458页 |
| 由GOCE引力梯度张量不变量确定卫星重力模型的半解析法;李建成 等;《武汉大学学报·信息科学版》;20160131;第41卷(第1期);第21-25页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107065025A (en) | 2017-08-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107065025B (en) | A kind of orbital elements estimation method based on gravimetric field gradient invariant | |
| CN109556632B (en) | INS/GNSS/polarization/geomagnetic integrated navigation alignment method based on Kalman filtering | |
| Chatfield | Fundamentals of high accuracy inertial navigation | |
| CN101344391B (en) | Lunar vehicle posture self-confirming method based on full-function sun-compass | |
| Crassidis et al. | New algorithm for attitude determination using Global Positioning System signals | |
| CN102175241B (en) | Autonomous astronomical navigation method of Mars probe in cruise section | |
| CN109556631B (en) | INS/GNSS/polarization/geomagnetic combined navigation system alignment method based on least squares | |
| CN112161632B (en) | Satellite formation initial positioning method based on relative position vector measurement | |
| CN101216319A (en) | Low orbit satellite multi-sensor fault tolerance autonomous navigation method based on federal UKF algorithm | |
| CN106871928A (en) | Strap-down inertial Initial Alignment Method based on Lie group filtering | |
| CN103674032A (en) | Satellite autonomous navigation system and method integrating pulsar radiation vector and timing observation | |
| CN111552003B (en) | Asteroid gravitational field full-autonomous measurement system and method based on ball satellite formation | |
| Liu et al. | X-ray pulsar/starlight Doppler integrated navigation for formation flight with ephemerides errors | |
| Wang et al. | Absolute navigation for Mars final approach using relative measurements of X-ray pulsars and Mars orbiter | |
| Xinlong et al. | An autonomous navigation scheme based on geomagnetic and starlight for small satellites | |
| Gu et al. | Optical/radio/pulsars integrated navigation for Mars orbiter | |
| Chen et al. | Gravity gradient tensor eigendecomposition for spacecraft positioning | |
| Kruger et al. | Observability analysis and optimization for angles-only navigation of distributed space systems | |
| Giorgi | Attitude determination | |
| CN116698048A (en) | Combined navigation method based on pulsar/inter-satellite ranging/landmark | |
| Chen et al. | Observability analysis for orbit determination using spaceborne gradiometer | |
| Pei et al. | Autonomous orbit determination using epoch-differenced gravity gradients and starlight refraction | |
| Sun et al. | On the feasibility of orbit determination from gravity gradient invariants | |
| Wang et al. | Autonomous orbit determination using pulsars and inter-satellite ranging for Mars orbiters | |
| Rosen | Analysis of hybrid satellite-to-satellite tracking and quantum gravity gradiometry architecture for time-variable gravity sensing missions |
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 |