CN106092105A - A kind of determination method of the strict regression orbit of near-earth satellite - Google Patents

Abstract

The determination method of the strict regression orbit of a kind of near-earth satellite, on the basis of rule of thumb formula obtains the orbital tracking discreet value of low order subgravity potential field situation Regression track, with semi-major axis of orbitaAnd orbit inclination angleiFor combining, according to semi-major axis of orbitaAnd orbit inclination angleiWith the relation of substar longitude and latitude, based on the Orbit simulation module of high-order subgravity potential field model, repeat to semi-major axis of orbitaAnd orbit inclination angleiIt is iterated revising, with eccentricityeAnd argument of perigeeωFor combination, for eccentricity vector properties of limit cycles, the method for average is used to repeat to eccentricityeAnd argument of perigeeωIt is iterated revising, until the regression accuracy of ascending node meets setting value.The present invention determines the strict regression orbit of near-earth satellite based on high-precision orbital dynamics, the track determining has higher regression accuracy for extraterrestrial target point, comparing traditional method based on low order subgravity potential field, high-precision dynamics of orbits is closer to reality, more using value.

Description

A kind of determination method of the strict regression orbit of near-earth satellite

Technical field

The invention belongs to spacecraft orbit dynamics technology field, particularly relate to the strict regression orbit of a kind of near-earth satellite Determine method

Background technology

After strict regression orbit requires one strict recursion period of experience, satellite can carry out high accuracy to extraterrestrial target point Revisit.For realizing the strict recurrence of track, the track product needed of design meets Sun synchronization repeating orbit and Frozen Orbit Characteristic.Wherein, it is optimized design according to Sun synchronization repeating orbit characteristic, it is possible to achieve revisiting of substar；Foundation is frozen Knot orbital characteristics is optimized design, it is possible to achieve the line of apsides stablizing in orbit plane, thus rail when ensureing that substar revisits The uniformity of road height.

Traditional regression orbit determines that method is based on low order subgravity potential field, and its major defect is that regression accuracy is not high, General at about 10km.

Content of the invention

The present invention provides the determination method of the strict regression orbit of a kind of near-earth satellite, comes really based on high-precision orbital dynamics Determining the strict regression orbit of near-earth satellite, the track of determination has higher regression accuracy for extraterrestrial target point, compares traditional Based on the method for low order subgravity potential field, high-precision dynamics of orbits is closer to reality, more using value.

In order to achieve the above object, the present invention provides the determination method of the strict regression orbit of a kind of near-earth satellite, comprise with Lower step: obtain the basis of the orbital tracking discreet value of low order subgravity potential field situation Regression track at rule of thumb formula On, with semi-major axis of orbit a and orbit inclination angle i for combination, according to semi-major axis of orbit a and orbit inclination angle i and substar longitude and latitude Relation, is derived by correction formula, and obtains iterative correction methods based on the Orbit simulation module of high-order subgravity potential field model, Repeat to be iterated semi-major axis of orbit a and orbit inclination angle i to revise, with eccentric ratio e and argument of perigee ω for combination, for partially The properties of limit cycles that the dynamic system of heart rate vector is had, uses the method for average to repeat to eccentric ratio e and argument of perigee ω It is iterated revising, it is achieved the freezing characteristic of track, until the regression accuracy of ascending node meets setting value.

Described regression accuracy is better than 5m.

Described rule of thumb formula obtains the orbital tracking discreet value a of low order subgravity potential field situation Regression track₀, i₀, e₀, ω₀Step comprise:

The given strict cycle T returning and corresponding track number of turns N, the orbital period of every railOnly consider low order Subgravity potential field situation, the discreet value of semi-major axis of orbit a is:

$a_{J_1}={\left(\frac{P}{2\mathrm{\π}}\sqrt{{\mathrm{GM}}_{\⊕}}\right)}^{\frac{2}{3}}$

$a_{J_2}=a_{J_1}+\frac{1}{J_2{\mathrm{GM}}_{\⊕}}{\left(\frac{4\stackrel{\·}{\mathrm{\Ω}}{a}_{J_1}^3}{3R_{\⊕}}\right)}^2-\frac{J_2{R}_{\⊕}^2}{a_{J_1}};$

Wherein,For Gravitational coefficient of the Earth,For earth radius；Subscript J of semi-major axis discreet value₁Represent that track moves Mechanics only considers disome situation, subscript J₂Represent and consider J₂Item terrestrial gravitation potential field；

The rate of change of right ascension of ascending node Ω meets:

$\frac{2\mathrm{\π}}{365.24\×86400}=\stackrel{\·}{\mathrm{\Ω}}=-\frac{3}{2}\sqrt{{\mathrm{GM}}_{\⊕}}{J}_2\frac{R_{\⊕}^2}{a_{J_2}^{3.5}}\cos i_{J_2};$

The discreet value of orbit inclination angle i is:

$i_{J_2}=\arccos (-\frac{2}{3}\frac{\stackrel{\·}{\mathrm{\Ω}}{a}_{J_2}^{3.5}}{\sqrt{{\mathrm{GM}}_{\⊕}}{J}_2{R}_{\⊕}^2});$

According to the requirement of Frozen Orbit, eccentric ratio e and argument of perigee ω meet:

The described step being iterated to semi-major axis of orbit a and orbit inclination angle i revising specifically comprises the steps of

The correction formula of step S2.1, derivation semi-major axis of orbit a and orbit inclination angle i

Step S2.2, the iterated revision formula obtaining semi-major axis of orbit a and orbit inclination angle i；

Assume longitude and latitude and orbital tracking meet functional relation λ=f (a, i),Obtain semi-major axis of orbit a and The iterated revision formula of orbit inclination angle i:

Step S2.3, according to semi-major axis of orbit discreet value a₀With orbit inclination angle discreet value i₀Calculate semi-major axis of orbit initial wink RadicalWith orbit inclination angle initial wink radical

Step S2.4, ascending node position determination module are according to semi-major axis initial wink radicalWith orbit inclination angle initial wink radicalCalculated the initial position r of ascending node by iterative approach₀With velocity v₀；

Step S2.5, the initial position r according to ascending node for the Orbit simulation module using high-order subgravity potential field model₀With Velocity v₀Carry out Orbit simulation, obtain the longitude and latitude difference Δ λ being separated by between two ascending nodes of a strict recursion period,

Step S2.6, by the longitude and latitude difference Δ λ between two ascending nodes,Substitute into semi-major axis of orbit a and orbit inclination angle i Iterated revision formula, obtain semi-major axis of orbit correction value Δ a and orbit inclination angle correction value Δ i；

Step S2.7, the longitude and latitude difference Δ λ judging between two ascending nodes,And semi-major axis of orbit correction value Δ a and Whether orbit inclination angle correction value Δ i meets threshold value simultaneously, if it is satisfied, by current semi-major axis of orbit correction value Δ a and track Inclination correction value Δ i, as final correction value, if be unsatisfactory for, carries out step S2.8；

Step S2.8, the semi-major axis of orbit correction value Δ a obtaining step S2.6 and orbit inclination angle correction value Δ i are as repeatedly For revised semi-major axis of orbit initial mean element a₀With orbit inclination angle initial mean element i₀, calculate the track after iterated revision half Major axis wink radicalWith orbit inclination angle wink radicalCarry out step S2.4.

The step of the correction formula of described derivation semi-major axis of orbit a and orbit inclination angle i specifically comprises:

Substar longitude and latitude meets:

Wherein ω_e=7.2921158 × 10^-5Rad/s, S₀Sidereal time for initial time Greenwich；

$\mathrm{\Ω}\≈{\mathrm{\Ω}}_0+\stackrel{\·}{\mathrm{\Ω}}(t-t_0)={\mathrm{\Ω}}_0-\frac{3J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{2{(1-e^2)}^2}\frac{\cos i}{a^{7/2}}(t-t_0);$

The finite term progression approximation of ascending node argument u meets:

$u\≈\mathrm{\ω}+M+(2e-\frac{e^3}{4})\sin M+\frac{5}{4}{e}^2\sin 2M+\frac{13}{12}{e}^{3}sin3M$

Wherein

Ascending node argument u with regard to the partial derivative of semi-major axis a is:

$\begin{array}{c}\frac{\∂u}{\∂a}=(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)\frac{\∂M}{\∂a}\\ =-\frac{3}{2}\sqrt{\frac{\mathrm{\μ}}{a^5}}(t-t_0)(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)\end{array};$

Ask f (a, i), g (a, i) with regard to a, the partial derivative of i, obtain:

$\left\{\begin{array}{c}\frac{\∂f}{\∂a}=\frac{21J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{4{(1-e^2)}^2}\frac{\cos i}{a^{9/2}}(t-t_0)+\frac{1}{1+{\left(\tan u\cos i\right)}^2}\frac{\cos i}{{\left(\cos u\right)}^2}\frac{\∂u}{\∂a}\\ \frac{\∂f}{\∂i}=-\frac{3J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{2{(1-e^2)}^2}\frac{\sin i}{a^{7/2}}(t-t_0)-\frac{1}{1+{\left(\tan u\cos i\right)}^2}\tan u\sin i\\ \frac{\∂g}{\∂a}=\frac{\cos u\sin i}{\sqrt{1-{\left(\sin u\sin i\right)}^2}}\frac{\∂u}{\∂a}\\ \frac{\∂g}{\∂i}=\frac{1}{\sqrt{1-{\left(\sin u\sin i\right)}^2}}\sin u\cos i\end{array}\right.;$

Value u=0 at ascending node, (t-t₀) value strict recursion period T, the correction formula of orbital tracking can be reduced to

$\left\{\begin{array}{c}\frac{\∂f}{\∂a}=\frac{21J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{4{(1-e^2)}^2}\frac{\cos i}{a^{9/2}}T-\frac{3}{2}\frac{\sqrt{\mathrm{\μ}}\cos i}{a^{5/2}}(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)T\\ \frac{\∂f}{\∂i}=-\frac{3J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{2{(1-e^2)}^2}\frac{\sin i}{a^{7/2}}T\\ \frac{\∂g}{\∂a}=-\frac{3}{2}\frac{\sqrt{\mathrm{\μ}}\sin i}{a^{5/2}}(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)T\\ \frac{\∂g}{\∂i}=0\end{array}\right..$

The step being iterated to eccentric ratio e and argument of perigee ω in described step S4 revising specifically comprises following step Rapid:

Step S4.1, the eccentricity vector e gathering multiple strict recursion period_x=e cos ω, e_y=e sin ω；

Step S4.2, the initial mean element as iteration next time for the average adding up eccentricity vector；

Whether the deviation between the average of the eccentricity vector obtaining twice before and after step S4.3, judgement is less than threshold value, as Fruit is, using current eccentricity correction value and argument of perigee correction value as final correction value, if it does not, carry out step S4.1。

Gather the eccentricity vector e of the strict recursion period of 4 months_x=e cos ω, e_y=e sin ω.

The step of the initial mean element as iteration next time for the average of described statistics eccentricity vector specifically comprises: utilize The eccentricity vector e collecting_x, e_yMapping, makes eccentricity vector form an approximation in the track Guan Bi of its variable space " justify ", using current " center of circle " as the initial value of the eccentric ratio e of next iteration and argument of perigee ω.

The present invention determines the strict regression orbit of near-earth satellite based on high-precision orbital dynamics, and the track of determination is for sky Between impact point there is higher regression accuracy, compare traditional method based on low order subgravity potential field, high-precision track move Mechanics is closer to reality, more using value.

Brief description

Fig. 1 is the flow chart of the determination method of the strict regression orbit of a kind of near-earth satellite that the present invention provides.

Fig. 2 is the iterative correction methods flow chart of the semi-major axis of orbit that provides of the present invention and orbit inclination angle.

Fig. 3 is the iterative correction methods flow chart of the eccentricity vector based on STK statistics that the present invention provides.

Fig. 4 is the design sketch of the iterated revision process of the eccentricity vector that the present invention provides.

Detailed description of the invention

Below according to Fig. 1～Fig. 4, illustrate presently preferred embodiments of the present invention.

As it is shown in figure 1, the present invention provides the determination method of the strict regression orbit of a kind of near-earth satellite, comprise the steps of

Step S1, rule of thumb formula obtain the orbital tracking discreet value of low order subgravity potential field situation Regression track (comprising semi-major axis of orbit a, orbit inclination angle i, eccentric ratio e and argument of perigee ω)；

Step S2, semi-major axis of orbit a and orbit inclination angle i are iterated revise；

Step S3, judge whether the regression accuracy of ascending node meets setting value, if so, then determine strict recurrence rail Road, if it is not, then carry out step S4；

In the present embodiment, regression accuracy is better than 5m；

Step S4, eccentric ratio e and argument of perigee ω are iterated revise, carry out step S2.

In described step S1, rule of thumb formula obtains the orbital tracking of low order subgravity potential field situation Regression track Discreet value a₀, i₀, e₀, ω₀Step comprise:

The given strict cycle T returning and corresponding track number of turns N, the orbital period of every railIf only considering low Order gravity potential field situation, the discreet value of semi-major axis of orbit a is:

$a_{J_1}={\left(\frac{P}{2\mathrm{\π}}\sqrt{{\mathrm{GM}}_{\⊕}}\right)}^{\frac{2}{3}}$

$a_{J_2}=a_{J_1}+\frac{1}{J_2{\mathrm{GM}}_{\⊕}}{\left(\frac{4\stackrel{\·}{\mathrm{\Ω}}{a}_{J_1}^3}{3R_{\⊕}}\right)}^2-\frac{J_2{R}_{\⊕}^2}{a_{J_1}};$

Wherein,For Gravitational coefficient of the Earth,For earth radius；Subscript J of semi-major axis discreet value₁Represent track power Learn and only consider disome situation, subscript J₂Represent and consider J₂Item terrestrial gravitation potential field；

The rate of change of right ascension of ascending node Ω meets:

$\frac{2\mathrm{\π}}{365.24\×86400}=\stackrel{\·}{\mathrm{\Ω}}=-\frac{3}{2}\sqrt{{\mathrm{GM}}_{\⊕}}{J}_2\frac{R_{\⊕}^2}{a_{J_2}^{3.5}}\cos i_{J_2};$

The discreet value of orbit inclination angle i is:

$i_{J_2}=\arccos (-\frac{2}{3}\frac{\stackrel{\·}{\mathrm{\Ω}}{a}_{J_2}^{3.5}}{\sqrt{{\mathrm{GM}}_{\⊕}}{J}_2{R}_{\⊕}^2});$

According to the requirement of Frozen Orbit, eccentric ratio e and argument of perigee ω meet:

As in figure 2 it is shown, in described step S2, be iterated the step tool revised to semi-major axis of orbit a and orbit inclination angle i Body comprises the steps of

The correction formula of step S2.1, derivation semi-major axis of orbit a and orbit inclination angle i

Substar longitude and latitude meets:

Wherein earth spin angle velocity ω_e=7.2921158 × 10^-5Rad/s, S₀Fixed star for initial time Greenwich When；

$\mathrm{\Ω}\≈{\mathrm{\Ω}}_0+\stackrel{\·}{\mathrm{\Ω}}(t-t_0)={\mathrm{\Ω}}_0-\frac{3J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{2{(1-e^2)}^2}\frac{\cos i}{a^{7/2}}(t-t_0);$

The finite term progression approximation of ascending node argument u meets:

$u\≈\mathrm{\ω}+M+(2e-\frac{e^3}{4})\sin M+\frac{5}{4}{e}^{2}sin2M+\frac{13}{12}{e}^{3}sin3M$

Wherein

Ascending node argument u with regard to the partial derivative of semi-major axis a is:

$\begin{array}{c}\frac{\∂u}{\∂a}=(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)\frac{\∂M}{\∂a}\\ =-\frac{3}{2}\sqrt{\frac{\mathrm{\μ}}{a^5}}(t-t_0)(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)\end{array};$

Ask f (a, i), g (a, i) with regard to a, the partial derivative of i, obtain:

$\left\{\begin{array}{c}\frac{\∂f}{\∂a}=\frac{21J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{4{(1-e^2)}^2}\frac{\cos i}{a^{9/2}}(t-t_0)+\frac{1}{1+{\left(\tan u\cos i\right)}^2}\frac{\cos i}{{\left(\cos u\right)}^2}\frac{\∂u}{\∂a}\\ \frac{\∂f}{\∂i}=-\frac{3J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{2{(1-e^2)}^2}\frac{\sin i}{a^{7/2}}(t-t_0)-\frac{1}{1+{\left(\tan u\cos i\right)}^2}\tan u\sin i\\ \frac{\∂g}{\∂a}=\frac{\cos u\sin i}{\sqrt{1-{\left(\sin u\sin i\right)}^2}}\frac{\∂u}{\∂a}\\ \frac{\∂g}{\∂i}=\frac{1}{\sqrt{1-{\left(\sin u\sin i\right)}^2}}\sin u\cos i\end{array}\right.;$

Value u=0 at ascending node, (t-t₀) value strict recursion period T, the correction formula of orbital tracking can be reduced to

$\left\{\begin{array}{c}\frac{\∂f}{\∂a}=\frac{21J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{4{(1-e^2)}^2}\frac{\cos i}{a^{9/2}}T-\frac{3}{2}\frac{\sqrt{\mathrm{\μ}}\cos i}{a^{5/2}}(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)T\\ \frac{\∂f}{\∂i}=-\frac{3J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{2{(1-e^2)}^2}\frac{\sin i}{a^{7/2}}T\\ \frac{\∂g}{\∂a}=-\frac{3}{2}\frac{\sqrt{\mathrm{\μ}}\sin i}{a^{5/2}}(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)T\\ \frac{\∂g}{\∂i}=0\end{array}\right.;$

Step S2.2, the iterated revision formula obtaining semi-major axis of orbit a and orbit inclination angle i；

Assume longitude and latitude and orbital tracking meet functional relation λ=f (a, i),Obtain semi-major axis of orbit a and The iterated revision formula of orbit inclination angle i:

Step S2.3, according to semi-major axis of orbit discreet value a₀(i.e. initial mean element) and orbit inclination angle discreet value i₀(i.e. initial Mean element) calculate semi-major axis of orbit initial wink radicalWith orbit inclination angle initial wink radical

Step S2.4, ascending node position determination module are according to semi-major axis initial wink radicalWith orbit inclination angle initial wink radicalCalculated the initial position r of ascending node by iterative approach₀With velocity v₀；

Step S2.5, the initial position r according to ascending node for the Orbit simulation module using high-order subgravity potential field model₀With Velocity v₀Carry out Orbit simulation, obtain the longitude and latitude difference Δ λ being separated by between two ascending nodes of a strict recursion period,

Step S2.6, by the longitude and latitude difference Δ λ between two ascending nodes,Substitute into semi-major axis of orbit a and orbit inclination angle i Iterated revision formula, obtain semi-major axis of orbit correction value Δ a and orbit inclination angle correction value Δ i；

Step S2.7, the longitude and latitude difference Δ λ judging between two ascending nodes,And semi-major axis of orbit correction value Δ a and Whether orbit inclination angle correction value Δ i meets such as lower threshold value simultaneously, if it is satisfied, by current semi-major axis of orbit correction value Δ a with Orbit inclination angle correction value Δ i, as final correction value, if be unsatisfactory for, carries out step S2.8；

||Δa||≤ε_a, or | | Δ i | |≤ε_i(||Δλ||≤ε_λ, or)；

Wherein, ε_aTake 0.05m, ε_iTake 0.001 °, ε_λTake (1.5 × 10^-6) °,Take (1.5 × 10^-6)°；

Step S2.8, the semi-major axis of orbit correction value Δ a obtaining step S2.6 and orbit inclination angle correction value Δ i are as repeatedly For revised semi-major axis of orbit initial mean element a₀With orbit inclination angle initial mean element i₀, calculate the track after iterated revision half Major axis wink radicalWith orbit inclination angle wink radicalCarry out step S2.4.

As it is shown on figure 3, the step being iterated to eccentric ratio e and argument of perigee ω in described step S4 revising is concrete Comprise the steps of

Step S4.1, the eccentricity vector e gathering multiple strict recursion period_x=e cos ω, e_y=e sin ω；

The dynamic system of eccentricity vector has " chummage limit cycle ", and for sun-synchronous orbit, eccentricity vector exists The period of change of its variable space is about 4 months, therefore, gathers the eccentricity vector e of the strict recursion period of 4 months_x=e Cos ω, e_y=e sin ω；

Step S4.2, the initial mean element as iteration next time for the average adding up eccentricity vector；

As shown in Figure 4, the eccentricity vector e collecting is utilized_x, e_yMapping, makes eccentricity vector at the rail of its variable space Mark Guan Bi forms " justifying " of an approximation, and Frozen Orbit requires eccentricity vector e_x, e_yAmplitude of variation as far as possible little, i.e. " justify " " radius " is little as far as possible, therefore using current " center of circle " (average) as the eccentric ratio e of next iteration and argument of perigee The initial value of ω；

Whether the deviation between the average of the eccentricity vector obtaining twice before and after step S4.3, judgement is less than threshold value (this In embodiment, threshold value is 10^-5), freeze if it is, illustrate that current eccentricity correction value and argument of perigee correction value meet Current eccentricity correction value and argument of perigee correction value are made by the freezing characteristic (i.e. eccentricity vector keeps constant) of track For final correction value, if it does not, carry out step S4.1.

In the present embodiment, the strict recursion period that design input is track 7 days, corresponding 101 orbital periods.Rule of thumb Formula, the initial estimate of available orbital tracking as shown in table 1.Described Orbit simulation module uses based on Matlab's Orbit simulation module, the 90*90 order gravitational potential field model choosing EGM2008 carries out Orbit simulation, and the initial of Orbit simulation is gone through 0 point 0 second when unit is 1 day 0 October in 2015, the initial simulation step length of Orbit simulation takes 5 seconds, every time the contracting of encryption gathering simulation step-length Being kept to previous 1/100, the position of ascending node determines encrypts collection twice, and the simulation step length that last encryption gathers is 5.0×10^-4Second.Use the STK data report functional realiey of STK software to the eccentric ratio e of multiple strict recursion periods and near-earth The collection of some argument ω, the coordinate system selection J2000 inertial coodinate system of Orbit simulation and STK data acquisition, kinetic model is only Consider terrestrial gravitation potential field.

As shown in table 1, initial estimate is exactly the orbital tracking discreet value obtaining in step S1, to track in step S2.7 The combination of semi-major axis and orbit inclination angle obtains Sun synchronization repeating orbit after being iterated revising, in step S4.3 to eccentricity and The combination of argument of perigee obtains Frozen Orbit after being iterated revising, and repeats to combine semi-major axis of orbit a and orbit inclination angle i, Eccentric ratio e and argument of perigee ω combination are iterated revising, until regression accuracy meets design and requires, obtain one group and strictly return Return orbit parameter.

The track mean element (0 point 0 second during initial 1 day 0 October of 2015 epoch) of table 1 each link correction gained

The present invention determines the strict regression orbit of near-earth satellite based on high-precision orbital dynamics, and the track of determination is for sky Between impact point there is higher regression accuracy, compare traditional method based on low order subgravity potential field, high-precision track move Mechanics is closer to reality, more using value.Although present disclosure has made detailed Jie by above preferred embodiment Continue, but it should be appreciated that the description above is not considered as limitation of the present invention.On those skilled in the art have read After stating content, multiple modification and replacement for the present invention all will be apparent from.Therefore, protection scope of the present invention should be by Appended claim limits.

Claims (8)

1. the determination method of the strict regression orbit of a near-earth satellite, it is characterised in that comprise the steps of rule of thumb public On the basis of formula obtains the orbital tracking discreet value of low order subgravity potential field situation Regression track, with semi-major axis of orbit a and rail Road inclination i is combination, according to the relation of semi-major axis of orbit a and orbit inclination angle i and substar longitude and latitude, is derived by revising public affairs Formula, and obtain iterative correction methods based on the Orbit simulation module of high-order subgravity potential field model, repeat to semi-major axis of orbit a and Orbit inclination angle i is iterated revising, with eccentric ratio e and argument of perigee ω for combination, for the dynamics system of eccentricity vector Unite had properties of limit cycles, use the method for average to repeat to be iterated to eccentric ratio e and argument of perigee ω revising, it is achieved rail The freezing characteristic in road, until the regression accuracy of ascending node meets setting value.

2. the determination method of the strict regression orbit of near-earth satellite as claimed in claim 1, it is characterised in that described recurrence essence Degree is less than 5m.

3. the determination method of the strict regression orbit of near-earth satellite as claimed in claim 2, it is characterised in that described according to warp Test formula and obtain the orbital tracking discreet value a of low order subgravity potential field situation Regression track₀, i₀, e₀, ω₀Step comprise:

The given strict cycle T returning and corresponding track number of turns N, the orbital period of every railOnly consider low order subgravity Potential field situation, the discreet value of semi-major axis of orbit a is:

$a_{J_1}={\left(\frac{P}{2\mathrm{\π}}\sqrt{{\mathrm{GM}}_{\⊕}}\right)}^{\frac{2}{3}}$

$a_{J_2}=a_{J_1}+\frac{1}{J_2{\mathrm{GM}}_{\⊕}}{\left(\frac{4\stackrel{\·}{\mathrm{\Ω}}{a}_{J_1}^3}{3R_{\⊕}}\right)}^2-\frac{J_2{R}_{\⊕}^2}{a_{J_1}};$

Wherein,For Gravitational coefficient of the Earth,For earth radius；Subscript J of semi-major axis discreet value₁Represent dynamics of orbits only Consider disome situation, subscript J₂Represent and consider J₂Item terrestrial gravitation potential field；

The rate of change of right ascension of ascending node Ω meets:

$\frac{2\mathrm{\π}}{365.24\×86400}=\stackrel{\·}{\mathrm{\Ω}}=-\frac{3}{2}\sqrt{{\mathrm{GM}}_{\⊕}}{J}_2\frac{R_{\⊕}^2}{a_{J_2}^{3.5}}\cos i_{J_2};$

The discreet value of orbit inclination angle i is:

$i_{J_2}=arccos(-\frac{2}{3}\frac{\stackrel{\·}{\mathrm{\Ω}}{a}_{J_2}^{3.5}}{\sqrt{{\mathrm{GM}}_{\⊕}}{J}_2{R}_{\⊕}^2});$

According to the requirement of Frozen Orbit, eccentric ratio e and argument of perigee ω meet:

ω=90 °.

4. the determination method of the strict regression orbit of near-earth satellite as claimed in claim 3, it is characterised in that described to track The step that semi-major axis a and orbit inclination angle i are iterated revising specifically comprises the steps of

The correction formula of step S2.1, derivation semi-major axis of orbit a and orbit inclination angle i

Step S2.2, the iterated revision formula obtaining semi-major axis of orbit a and orbit inclination angle i；

Assume longitude and latitude and orbital tracking meet functional relation λ=f (a, i),Obtain semi-major axis of orbit a and track The iterated revision formula of inclination angle i:

Step S2.3, according to semi-major axis of orbit discreet value a₀With orbit inclination angle discreet value i₀Calculate semi-major axis of orbit initial wink radicalWith orbit inclination angle initial wink radical

Step S2.4, ascending node position determination module are according to semi-major axis initial wink radicalWith orbit inclination angle initial wink radicalLogical Cross iterative approach and calculate the initial position r of ascending node₀With velocity v₀；

Step S2.5, the initial position r according to ascending node for the Orbit simulation module using high-order subgravity potential field model₀And speed Vector v₀Carry out Orbit simulation, obtain the longitude and latitude difference Δ λ being separated by between two ascending nodes of a strict recursion period,

Step S2.6, by the longitude and latitude difference Δ λ between two ascending nodes,Substitute into semi-major axis of orbit a and orbit inclination angle i repeatedly For correction formula, obtain semi-major axis of orbit correction value Δ a and orbit inclination angle correction value Δ i；

Step S2.7, the longitude and latitude difference Δ λ judging between two ascending nodes,And semi-major axis of orbit correction value Δ a and track Whether inclination correction value Δ i meets threshold value simultaneously, if it is satisfied, by current semi-major axis of orbit correction value Δ a and orbit inclination angle Correction value Δ i is as final correction value, if be unsatisfactory for, carries out step S2.8；

Step S2.8, the semi-major axis of orbit correction value Δ a obtaining step S2.6 and orbit inclination angle correction value Δ i repair as iteration Semi-major axis of orbit initial mean element a after just₀With orbit inclination angle initial mean element i₀, calculate the semi-major axis of orbit after iterated revision Wink radicalWith orbit inclination angle wink radicalCarry out step S2.4.

5. the determination method of the strict regression orbit of near-earth satellite as claimed in claim 4, it is characterised in that described derivation rail The step of the correction formula of road semi-major axis a and orbit inclination angle i specifically comprises:

Substar longitude and latitude meets:

Wherein ω_e=7.2921158 × 10^-5Rad/s, S₀Sidereal time for initial time Greenwich；

$\mathrm{\Ω}\≈{\mathrm{\Ω}}_0+\stackrel{\·}{\mathrm{\Ω}}(t-t_0)={\mathrm{\Ω}}_0-\frac{3J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{2{(1-e^2)}^2}\frac{\cos i}{a^{7/2}}(t-t_0);$

The finite term progression approximation of ascending node argument u meets:

$u\≈\mathrm{\ω}+M+(2e-\frac{e^3}{4})\sin M+\frac{5}{4}{e}^{2}sin2M+\frac{13}{12}{e}^3\sin 3M$

Wherein

Ascending node argument u with regard to the partial derivative of semi-major axis a is:

$\begin{array}{c}\frac{\∂u}{\∂a}=(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)\frac{\∂M}{\∂a}\\ =-\frac{3}{2}\sqrt{\frac{\mathrm{\μ}}{a^5}}(t-t_0)(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)\end{array};$

Ask f (a, i), g (a, i) with regard to a, the partial derivative of i, obtain:

$\left\{\begin{array}{c}\frac{\∂f}{\∂a}=\frac{21J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{4{(1-e^2)}^2}\frac{\cos i}{a^{9/2}}(t-t_0)+\frac{1}{1+{\left(\tan u\cos i\right)}^2}\frac{\cos i}{{\left(\cos u\right)}^2}\frac{\∂u}{\∂a}\\ \frac{\∂f}{\∂i}=-\frac{3J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{2{(1-e^2)}^2}\frac{\sin i}{a^{7/2}}(t-t_0)-\frac{1}{1+{\left(\tan u\cos i\right)}^2}\tan u\sin i\\ \frac{\∂g}{\∂a}=\frac{\cos u\sin i}{\sqrt{1-{\left(\sin u\sin i\right)}^2}}\frac{\∂u}{\∂a}\\ \frac{\∂g}{\∂i}=\frac{1}{\sqrt{1-{\left(\sin u\sin i\right)}^2}}\sin u\cos i\end{array}\right.;$

Value u=0 at ascending node, (t-t₀) value strict recursion period T, the correction formula of orbital tracking can be reduced to

$\left\{\begin{array}{c}\frac{\∂f}{\∂a}=\frac{21J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{4{(1-e^2)}^2}\frac{\cos i}{a^{9/2}}T-\frac{3}{2}\frac{\sqrt{\mathrm{\μ}}\cos i}{a^{5/2}}(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)T\\ \frac{\∂f}{\∂i}=\frac{3J_2{R}_{\⊕}^2\sqrt{\mathrm{\μ}}}{2{(1-e^2)}^2}\frac{\sin i}{a^{7/2}}T\\ \frac{\∂g}{\∂a}=-\frac{3}{2}\frac{\sqrt{\mathrm{\μ}}\sin i}{a^{5/2}}(1+(2e-\frac{e^3}{4})\cos M+\frac{5}{2}{e}^2\cos 2M+\frac{13}{4}{e}^3\cos 3M)T\\ \frac{\∂g}{\∂i}=0\end{array}\right..$

6. the determination method of the strict regression orbit of near-earth satellite as claimed in claim 5, it is characterised in that described step S4 In the step that is iterated to eccentric ratio e and argument of perigee ω revising specifically comprise the steps of

Step S4.1, the eccentricity vector gathering multiple strict recursion period

Step S4.2, the initial mean element as iteration next time for the average adding up eccentricity vector；

Whether the deviation between the average of the eccentricity vector obtaining twice before and after step S4.3, judgement is less than threshold value, if it is, Using current eccentricity correction value and argument of perigee correction value as final correction value, if it does not, carry out step S4.1.

7. the determination method of the strict regression orbit of near-earth satellite as claimed in claim 6, it is characterised in that gather 4 months The eccentricity vector e of strict recursion period_x=e cos ω, e_y=e sin ω.

8. the determination method of the strict regression orbit of near-earth satellite as claimed in claim 6, it is characterised in that described statistics is inclined The step of the initial mean element as iteration next time for the average of heart rate vector specifically comprises: utilize the eccentricity vector collecting e_x, e_yMapping, makes eccentricity vector form " the justifying " approximating in the track Guan Bi of its variable space, with current " center of circle " Initial value as the eccentric ratio e of next iteration and argument of perigee ω.