CN103514362A - Two-line element generation method based on model error compensation - Google Patents

Abstract

The invention discloses a two-line element generation method based on model error compensation. According to the method, first, observed position and speed vectors are converted to a true equator mean equinox coordinate system through a geocentric equatorial inertial coordinate system, then corresponding instant orbital elements are computed and are used as the initial estimated values of two-line element iterative computation; the position and speed vectors at the corresponding sampling instant are respectively predicted according to the initial estimated values, the difference values between the position and speed vectors and the position and speed vectors at the practical observing sampling instant are computed, whether a termination condition is met is judged; through a numerical differentiation method, a partial derivative matrix about the two-line elements of a simplified pervasive perturbation orbit prediction function at the sampling instant is computed; by solving a compensation least square problem, the correction of the two-line elements is obtained; new two-line elements are computed, and repeated iterative operation is carried out until the accuracy requirement is met. According to the method, the situations of single-point fitting and sampling fitting can be processed, and average errors are small when the two-line elements obtained by fitting are used for orbit prediction.

Description

The two row radical generation methods based on Compensation for Model Errors

Technical field

The invention belongs to extraterrestrial target track computing technique, particularly a kind of two row radical generation methods based on Compensation for Model Errors.

Background technology

Two row orbital trackings (TLE) are the forms that the US Strategic Command issue extraterrestrial target catalogue data adopts, be a kind of mean orbit unit based on Kepler's elements, be widely used in describing the extraterrestrial target flight track of (comprising moonlet and space junk etc.).For two row radicals, it is with ad hoc approach, to have removed the median orbital elements of periodic perturbation item, it must make for forecasting position and the speed in a certain moment of extraterrestrial target together with specific orbit prediction model, for this forecasting model of near-earth target, is exactly the SGP4 model that North American Air Defense Command (NORAD) provides.As the class analytic model being most widely used at present, it has considered the impact of the perturbation factors such as the non-spherical J2 of the earth, J3, the perturbation of J4 item, atmospheric perturbation and lunisolar gravitational perturbation, has made good balance between track forecast precision and counting yield.

At present, two row radicals of most of extraterrestrial target can regularly be announced renewal in Goddard space research center, as long as by TLE substitution disclosed SGP4 model, can obtain position and the speed of any time extraterrestrial target, but TLE generating algorithm is still unexposed so far, this makes the further application based on TLE be subject to larger restriction.Walter(1.Walter H G.Conversion of osculating obbital elements into mean elements[J as far back as Germany in 1967] .The Astronomical Journal, 1967,72 (8): 994-997) just under the orbit theory meaning of Brouwer and Kozai, studied by instantaneous orbit radical and be converted to corresponding median orbital elements problem, proposed an iterative algorithm and proved its convergence.After SGP4 model occurs, as a kind of median orbital elements special and that be widely used, the generating algorithm of two row radicals has caused research widely.Wherein, (the 2.Jochim E F such as Jochim, Gill E, Montenbruck O, et al.GPS based on-board and onground orbit operations for small satellites[J] .Acta Astronautica, 1996,39 (9): 917-922) usining the Onboard GPS Data of extraterrestrial target generates the TLE of this target in real time as the defeated people of observation, but the method is only applicable to cooperative target.Lee(3.Byoung S L.NORAD TLE conversion from osculating orbital elements[J] .Journal of Astronomy and Space Sciences, 2002,19 (4): 395-402) proposition utilizes the instantaneous osculating orbit data fitting of extraterrestrial target to go out the TLE of this target, but only studied the method for single-point matching, and do not considered the impact of the atmospherical drag factor; Dwight(4.Dwight E A.Computing NORAD mean orbital elements from a state vector[D] .Ohio:Air Force Institute of Technology, Wright-Patterson Air Force Base, 1994) consider that Liao You position and velocity generate the problem of two row radicals, Newton method and direct two kinds of methods of iteration have been proposed, but the method unresolvedly occur that singular point causes the problem that cannot restrain in excentricity and orbit inclination when very little; (5. Liu is bright for Liu Guangming etc., Wen Yuanlan, Liao Ying. the duplicate rows orbital tracking generating algorithm [J] based on without singular transformation. systems engineering and electronic technology, 2011,33 (5): 1104-1107) proposed the extraterrestrial target TLE generating algorithm based on without singular transformation, for the problem that may occur singular point in TLE sampling fit procedure, introduce without singular orbit radical and replace Kepler's elements and provided the approximate analysis expression formula of dbjective state vector about two row radical partial derivative matrixs, but the objective function of its matching has only been considered position vector.(the 6. Hu Min such as Hu Min, CENG GUOQIANG. the exchange of median orbital elements and osculating element [J]. Spacecraft TT&C journal, 2012,31 (2): 77-81) based on TLE and SGP4 model, from the pros and cons, summed up the interchange issues of median orbital elements and osculating element, only with single-point fitting process, calculated.

Because SGP4 model itself has error, this has affected the precision of two row radical matchings, the problem existing for current research, the present invention has proposed the error compensation Optimized model of two row radical matchings on the basis of considering SGP4 model error, by unified processing of situation of single-point matching and sampling matching, while making two row radicals that matching obtains for orbit prediction, average error is less.

Summary of the invention

Goal of the invention: the object of the present invention is to provide a kind of extraterrestrial target two row radical generation methods based on simplifying pervasive perturbed orbit forecasting model (SGP4).

Technical scheme: the technical solution that realizes the object of the invention is: a kind of two row radical generation methods based on Compensation for Model Errors, step is as follows:

Step 1, the position observing and velocity are converted into true equator mean equinox (TEME) coordinate system by equator, the earth's core inertia (J2000) coordinate system;

Position and the velocity of step 2, utilization sampling initial time are calculated corresponding instantaneous orbit radical the initial estimate using it as two row radical iterative computation;

Step 3, utilization are simplified pervasive perturbed orbit forecasting model SGP4(and are got B ^*=0) initial estimate of and two row radicals forecasts respectively position and the velocity vector of corresponding sampling instant, calculates itself and the position of actual observation sampling instant and the difference of velocity vector and judges whether to meet end condition.If meet, algorithm stops; If do not meet, go to step 4;

Step 4, by the pervasive perturbed orbit predictor of numerical differentiation computational short cut in corresponding sampling instant the partial derivative matrix about two row radicals;

Step 5, choose suitable smoothing factor and regular matrix, by solving compensation least square problem, obtain the correction of two row radicals;

Step 6, calculate two new row radicals, and carry out iteration, until meet accuracy requirement.

Beneficial effect: compared with prior art, its remarkable advantage is in the present invention: 1) simplify on the model error basis of pervasive perturbed orbit forecasting model having considered, can unify to process the situation of single-point matching and sampling matching; 2) objective function of matching has been taken into account position and velocity vector simultaneously, and while making two row radicals that matching obtains for orbit prediction, average error is less.

Accompanying drawing explanation

Fig. 1 is that the two row radicals based on Compensation for Model Errors generate method flow diagram.

Embodiment

Below in conjunction with Figure of description, the present invention is described in further detail:

As shown in Figure 1, the present invention relates to a kind of two row radical generation methods based on Compensation for Model Errors, concrete steps are as follows:

Step 1, by the position observing and velocity Y ₀, Y ₁..., Y _nby equator, the earth's core inertia (J2000) coordinate system, be converted into true equator mean equinox (TEME) coordinate system, do as down conversion:

${\left[\begin{array}{c}x\\ y\\ z\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\\ \stackrel{\·}{z}\end{array}\right]}_{\mathrm{TEME}}=R_z^{-1}\left(\mathrm{\Δ\φ}\cos \stackrel{\‾}{\mathrm{\ϵ}}\right)\·N\·P\·{\left[\begin{array}{c}x\\ y\\ z\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\\ \stackrel{\·}{z}\end{array}\right]}_{J2000},$

In formula,

represent position and velocity and Y that t observes constantly _i=Y (t _i), i=0,1 ..., n, R _zfor rotation of coordinate matrix, N is nutating matrix, and P is precession of the equinoxes matrix, and Δ φ is gold nutating,

for mean obliquity.

Step 2, utilization sampling initial time t ₀shi position r=(x, y, z) ^tand speed

calculate corresponding instantaneous orbit radical, the initial estimate X using it as two row radical iterative computation ₀=(e ₀, i ₀, Ω ₀, ω ₀, M ₀, n ₀) ^t; Be specially:

Step 2-1, to calculate intermediate variable as follows:

W _x=h _x/||h||，W _y=h _y/||h||,W _x=h _x/||h||，

Wherein, $G{M}_{\⊕}=398600.4415k{m}^3/s^2$ And

$h=\left(\begin{array}{c}{h}_x\\ h_y\\ h_z\end{array}\right)=\left(\begin{array}{c}y\stackrel{\·}{z}-z\stackrel{\·}{y}\\ z\stackrel{\·}{x}-x\stackrel{\·}{z}\\ x\stackrel{\·}{y}-y\stackrel{\·}{x}\end{array}\right),$ $\left|\right|h\left|\right|=\sqrt{h_x^2+h_y^2+h_z^2},$ $\left|\right|r\left|\right|=\sqrt{x^2+y^2+z^2},$ $\left|\right|\stackrel{\·}{r}\left|\right|=\sqrt{{\stackrel{\·}{x}}^2+{\stackrel{\·}{y}}^2+{\stackrel{\·}{z}}^2};$

Step 2-2, calculating orbit inclination and right ascension of ascending node are as follows:

$i_0=\arctan \left(\frac{\sqrt{{W_x}^2+{W_y}^2}}{W_z}\right),$ ${\mathrm{\Ω}}_0=\arctan (-\frac{W_x}{W_z});$

Step 2-3, calculating mean angular velocity and excentricity are as follows:

$n_0=\sqrt{\frac{G{M}_{\⊕}}{a^3}},$ $e_0=\sqrt{1-\frac{p}{a}};$

Step 2-4, to calculate mean anomaly as follows:

M ₀=E ₀-esin(E ₀),

Wherein,

$E_0=\arctan \left(\frac{(x\stackrel{\·}{x}+y\stackrel{\·}{y}+z\stackrel{\·}{z})/\left(a^2{n}_0\right)}{1-\left|\right|r\left|\right|/a}\right);$

Step 2-5, to calculate argument of perigee as follows:

${\mathrm{\ω}}_0=\arctan \left(\frac{z}{-x{W}_y+y{W}_x}\right)-\arctan \left(\frac{\sqrt{1-e_0^2}\sin \left(E_0\right)}{\cos \left(E_0\right)-e_0}\right).$

Step 3, utilization are simplified pervasive perturbed orbit forecasting model SGP4(and are got B ^*=0) the current estimated value of and two row radicals is forecast respectively position and the velocity vector of corresponding sampling instant, and calculates itself and the position of actual observation sampling instant and the difference of velocity vector; Be specially:

Step 3-1, utilize the predictor F (X of SGP4 model ₀, t ₀, t) calculate t ₀, t ₁..., t _nposition and velocity vector constantly

?

$\tilde{y}\left(t\right)=F(X_0,t_0,t),$ t=t ₀,t ₁,…,t _n.

Wherein, X ₀for t ₀the initial estimate of the moment two row radicals.

Step 3-2, calculate the position of current radical forecast and the position of velocity vector and actual observation and velocity vector at sampling instant t ₀, t ₁..., t _ndifference f ₀, f ₁..., f _n,

$f\left(t\right)=Y\left(t\right)-\tilde{Y}\left(t\right),$ t=t ₀,t ₁,…,t _n.

Step 3-3, judge whether difference size meets end condition.If meet, algorithm stops; If do not meet, go to step 4;

Step 4, by the pervasive perturbed orbit predictor of numerical differentiation computational short cut at sampling instant t ₀, t ₁..., t _npartial derivative matrix A about two row radicals _i=A (t _i), i=0,1 ..., n,

$A\left(t_i\right)=\frac{\∂F(X,t_0,t)}{\∂X}{|}_{t=t_i}={\left[\begin{array}{cccccc}\frac{\∂y_1}{\∂x_1}& \frac{\∂y_1}{\∂x_2}& \frac{\∂y_1}{\∂x_3}& \frac{\∂y_1}{\∂x_4}& \frac{\∂y_1}{\∂x_5}& \frac{\∂y_1}{\∂x_6}\\ \frac{\∂y_2}{\∂x_1}& \frac{\∂y_2}{\∂x_2}& \frac{\∂y_2}{\∂x_3}& \frac{\∂y_2}{\∂x_4}& \frac{\∂y_2}{\∂x_5}& \frac{\∂y_2}{\∂x_6}\\ \frac{\∂y_3}{\∂x_1}& \frac{\∂y_3}{\∂x_2}& \frac{\∂y_3}{\∂x_3}& \frac{\∂y_3}{\∂x_4}& \frac{\∂y_3}{\∂x_5}& \frac{\∂y_3}{\∂x_6}\\ \frac{\∂y_4}{\∂x_1}& \frac{\∂y_4}{\∂x_2}& \frac{\∂y_4}{\∂x_3}& \frac{\∂y_4}{\∂x_4}& \frac{\∂y_4}{\∂x_5}& \frac{\∂y_4}{\∂x_6}\\ \frac{\∂y_5}{\∂x_1}& \frac{\∂y_5}{\∂x_2}& \frac{\∂y_5}{\∂x_3}& \frac{\∂y_5}{\∂x_4}& \frac{\∂y_5}{\∂x_5}& \frac{\∂y_5}{\∂x_6}\\ \frac{\∂y_6}{\∂x_1}& \frac{\∂y_6}{\∂x_2}& \frac{\∂y_6}{\∂x_3}& \frac{\∂y_6}{\∂x_4}& \frac{\∂y_6}{\∂x_5}& \frac{\∂y_6}{\∂x_6}\end{array}\right]}_{t=t_i},i=\mathrm{0,1},\·\·\·,n.$

Wherein,

$F(X,t_0,t)=\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\\ y_5\\ y_6\end{array}\right],$ $X=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\\ x_6\end{array}\right],$

$\frac{\∂y_i}{\∂x_j}=\frac{y_i(x_1,\·\·\·,x_j+\mathrm{\Δ}{x}_j/2,\·\·\·,x_6,t)-y_i(x_1,\·\·\·,x_j-\mathrm{\Δ}{x}_j/2,\·\·\·,x_6,t)}{\mathrm{\Δ}{x}_j},i,j=\mathrm{1,2},\·\·\·,6.$

Δ x _jcan choose according to accuracy requirement general desirable Δ x _j=10 ^-3x _j.

Step 5, choose suitable smoothing factor and regular matrix, by solving compensation least square problem, obtain the correction of two row radicals; Be specially:

Step 5-1, choose smoothing factor α and regular matrix R=P ^tp, and compute matrix S=(I+ α R) ^-1, wherein

α can rule of thumb attempt choosing, and span is generally 0 < α < 1, the unit matrix that I is corresponding exponent number.

Step 5-2, calculate the correction vector of two row radicals:

ΔX ₀=(A ^T(I-S)A) ^-1A ^T(I-S)f.

Wherein,

$f=\left[\begin{array}{c}{f}_0\\ f_1\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ f_n\end{array}\right],$ $A=\left[\begin{array}{c}{A}_0\\ A_1\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ A_n\end{array}\right],$

Step 6, calculate two new row radicals, and carry out iteration, until meet accuracy requirement.

$X_0=X_0^{*}+\mathrm{\&Delta;}{X}_0,$ Can make $X_0^{*}=X_0$ Carry out iterative.

As from the foregoing, the present invention proposes about two row radicals and more generally optimize computation model, be both applicable to single-point matching and be also applicable to the matching of sampling, compare the existing method scope of application wider.Both considered the model error of SGP4, and also considered that the deviation of velocity was as optimization aim, and compared the two row radicals that existing method obtains, its mean accuracy for orbit prediction is higher.

Claims (7)

1. two row radical generation methods based on Compensation for Model Errors, is characterized in that, comprise the steps:

1) position observing and velocity are converted into true equator mean equinox coordinate system by equator, the earth's core inertial coordinates system;

2) position and the velocity of sampling initial time are converted to the initial estimate of instantaneous orbit radical using it as two row radical iterative computation;

3) initial estimate that utilize to simplify pervasive perturbed orbit forecasting model and two row radicals forecasts respectively position and the velocity vector of corresponding sampling instant, calculates itself and the position of actual observation sampling instant and the difference of velocity vector and judges whether to meet end condition.If meet, algorithm stops; If do not meet, go to step 4;

4) by the pervasive perturbed orbit predictor of numerical differentiation computational short cut in sampling instant the partial derivative matrix about two row radicals;

5) choose suitable smoothing factor and regular matrix, by solving compensation least square problem, obtain the correction of two row radicals;

6) calculate two new row radicals, and carry out iteration, until meet accuracy requirement.

2. the two row radical generation methods based on Compensation for Model Errors according to claim 1, is characterized in that, in step 1, to the specific formula for calculation of observation data coordinate transform, are:

In formula,

represent position and velocity and Y that t observes constantly _i=Y (t _i), i=0,1 ..., n., R _zfor rotation of coordinate matrix, N is nutating matrix, and P is precession of the equinoxes matrix, and △ φ is gold nutating,

for mean obliquity.

3. the two row radical generation methods based on Compensation for Model Errors according to claim 1, is characterized in that, step 2 is utilized sampling initial time t ₀shi position r=(x, y, z) ^tand speed

calculate corresponding instantaneous orbit radical, the initial estimate X using it as two row radical iterative computation ₀=(e ₀, i ₀, Ω ₀, ω ₀, M ₀, n ₀) ^t; Be specially:

1) calculating intermediate variable is as follows:

Wherein,

and

2) calculating orbit inclination and right ascension of ascending node are as follows:

3) calculating mean angular velocity and excentricity are as follows:

4) calculating mean anomaly is as follows:

M ₀=E ₀-esin(E ₀),

Wherein,

5) calculating argument of perigee is as follows:

4. the two row radical generation methods based on Compensation for Model Errors according to claim 1, is characterized in that, step 3 utilization is simplified pervasive perturbed orbit forecasting model SGP4(and got B ^*=0) the current estimated value of and two row radicals is forecast respectively position and the velocity vector of corresponding sampling instant, and calculates itself and the position of actual observation sampling instant and the difference of velocity vector, and concrete steps are:

1) utilize the predictor F (X of SGP4 model ₀, t ₀, t) calculate t ₀, t ₁..., t _nposition and velocity vector constantly

?

Wherein, X ₀for t ₀the initial estimate of the moment two row radicals.

2) calculate the position of current radical forecast and the position of velocity vector and actual observation and velocity vector at sampling instant t ₀, t ₁..., t _ndifference f ₀, f ₁..., f _n,

3) judge whether difference size meets end condition.If meet, algorithm stops; If do not meet, go to step 4.

5. the two row radical generation methods based on Compensation for Model Errors according to claim 1, is characterized in that, pass through the pervasive perturbed orbit predictor of numerical differentiation computational short cut at sampling instant t in step 4 ₀, t ₁..., t _npartial derivative matrix A about two row radicals _i=A (t _i), i=0,1 ..., n, specific formula for calculation is:

Wherein,

6. two row radical generation methods based on Compensation for Model Errors according to claim 1, it is characterized in that, in step 5, choose suitable smoothing factor and regular matrix, obtain the correction of two row radicals by solving compensation least square problem, concrete steps are:

1) choose smoothing factor α and regular matrix R=P ^tp, and compute matrix S=(I+ α R) ^-1, wherein

α can rule of thumb attempt choosing, and span is generally 0 < α < 1, the unit matrix that I is corresponding exponent number.

2) calculate the correction vector of two row radicals:

△X ₀=(A ^T(I-S)A) ^-1A ^T(I-S)f.

Wherein,

7. the two row radical generation methods based on Compensation for Model Errors according to claim 1, is characterized in that, calculate two new row radicals and carry out iteration and use following formula in step 6:

order

iterate and solve.

Images