Estimation of ballistic coefficients of space debris using the ratios between different objects

Zhejun LU,Weidong HU

College of Electronic Science and Engineering,National University of Defense Technology,Changsha 410073,China

Estimation of ballistic coefficients of space debris using the ratios between different objects

Zhejun LU*,Weidong HU

College of Electronic Science and Engineering,National University of Defense Technology,Changsha 410073,China

Available online 20 April 2017

*Corresponding author.

E-mail addresses:luzhejun@nudt.edu.cn(Z.LU),wdhu@nudt.edu.cn(W.HU).

Peer review under responsibility of Editorial Committee of CJA.

http://dx.doi.org/10.1016/j.cja.2017.03.009

1000-9361?2017 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.

This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

Production and hosting by Elsevier

This paper proposes a new method to estimate the ballistic coefficients(BC)of low earth orbit space debris.The data sources are the historical two-line elements(TLEs).Since the secular variation of semi-major axes is mainly caused by the drag perturbation for space objects with perigee altitude below 600 km,the ballistic coefficients are estimated based on variation of the mean semi-major axes derived from the TLEs.However,the approximate parameters used in the calculation have error,especially when the upper atmosphere densities are difficult to obtain and always estimated by empirical model.The proportional errors of the approximate parameters are cancelled out in the form of ratios,greatly mitigating the effects of model error.This method has been also been validated for space objects with perigee altitude higher than 600 km.The relative errors of estimated BC values from the new method are significantly smaller than those from the direct estimation methods used in numerical experiments.The estimated BC values are used for the prediction of the semi-major axes,and good performance is obtained.This process is also a feasible method for prediction over a long period of time without an orbital propagator model.

?2017 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

Ballistic coefficients;

Grouping;

Ratio;

Space debris;

TLE

1.Introduction

At operationally important orbits,there is a significantly increased amount of space debris created by spacecraft launch,loss,collision and explosion.1,2The presence of such debris causes a risk of collision,which threatens the safe operation of aircraft.3–5Management of the space debris population,which includes its cataloguing,prediction and mitigation,is crucial for the continued security of the space environment.For this reason,many methods have been proposed for space debris control and removal.6–9

Since space debris is consistently increasing,it is becoming ever more important to confirm its origin and assert clear legal responsibility.However,because of the large number and small size of space debris,cataloguing based on direct observation is difficult,and a ‘group catalogue‘is a more appropriate and efficient method.10The cataloguing of debris should be started immediately following its production,and debris derived from the same object should be processed together as a group.Space debris,more or less,maintains the orbit of its parent satellite.Thus,debris objects from the same parent have similar orbital elements and there is stable variable relationship between the values of their parameters.The ratios between these parameters provide new information for a debris group.Based on this,the state of individual debris pieces can be deduced from other debris in the same group.This information can be used as a criterion for judging which debris group a particular piece of debris belongs to,and the ‘group catalogue”can be achieved.

A ballistic coefficients(BC)is an important parameter for the research of space objects and consists of object’s drag coefficients(CD),mass and the cross sectional area in the direction of motion relative to the atmosphere.The atmospheric drag is a significant perturbing force for low earth orbit(LEO)space objects.The physical properties of the upper atmosphere and the ballistic coefficients of the space objects are needed for the calculation of the atmospheric force.The BC value is usually unknown for space debris and hard to measure.Thus,many methods have been proposed to estimate this parameter.

One way is to use the additional parameterB*given in the two-line elements sets(TLEs).11However,theB*is a fitting parameter in the process of producing the TLEs,and includes the biases related to model error.12Moreover,the value ofB*may even be negative in TLEs.13Thus,theB*is an inaccurate and unreliable estimation of the satellite’s ballistic coefficients.

Combined with the atmospheric drag equation,the BC value can be obtained from the filtering process.However,since stable and accurate measurements cannot be obtained,accurate target motion estimation is difficult to extract,and this process has very low accuracy.14

Using satellite track data through a 30-year historical time span,a batch least-squares differential correction algorithm is used to estimate the ballistic coefficients for the use of high accuracy satellite drag model(HASDM).15,16Since this method is,to a large extent,based on measurement data,it cannot be widely used for space objects.

The methods used in recent years have always estimated ballistic coefficients from long-term TLE data,e.g.the methods proposed by Picone et al.,17Saunders et al.18and Sang et al.19Based on the simplified relationship between mean motion and atmospheric drag,the atmospheric drag can be estimated through the mean motion extracted from the TLEs.The results are obtained with use of over 30 years of TLE data.However,since orbit propagator and empirical atmospheric models are used in those methods,accuracy depends on precise modeling and is limited.Consequently,of the quality of an estimated ballistic coefficients is dependent upon model accuracy.

From the works mentioned above,it can be found that the most difficult aspect of BC estimation is the availability of the accurate values for the parameters needed.Many cannot be measured and use empirical models for approximation.This paper presents a method to estimate the ratios of the ballistic coefficients between different objects.The proportional errors of the approximate parameters are cancelled out in the form of ratios,greatly mitigating the effects of the model error.Through this approach,the relative error of estimated BC values is significantly reduced compared with direct estimation methods.Moreover,it has been validated that this approach can be used for the objects with a perigee altitude higher than 600 km.The proposed method is used for space debris from a debris group,and results show stable ratios for estimated BC values between different debris objects.These ratios can be used to confirm that the debris objects originated from the same source and also to deduce the states of debris objects from the others in the same group.This process provides new information about the debris group and gives theoretical and methodological support for a group catalogue.The estimated BC values are also used for the prediction of semi-major axes without an orbital propagator model.

The paper is organized as follows:Section 2 proposes a comparison approach to estimate the ballistic coefficients.The fitting process for the change of semi-major axes is necessary and presented in Section 3.The results and analyses of numerical experiments are shown in Section 4.Section 5 tests the proposed method for the LEO objects higher than 600 km in perigee altitude.The method is used for space debris in Section 6.Section 7 discusses the inverse use of the proposed approach for orbital prediction.Finally,conclusions are drawn in Section 8.

2.Ballistic coefficients estimation

2.1.Atmospheric drag

Atmospheric drag is a nonconservation force and continuously affects orbit semi-major axes and orbit period decrease.20For LEO objects,atmospheric drag is the major source of the secular term on semi-major axes.The mean semi-major axis is a mean orbital element without periodic terms.21The variation of mean semi-major axes does not include periodic gravitational perturbation,and remaining secular gravitational terms are small.17Thus,the variation in meansemi-major axes mainly reflects the effects of atmospheric drag.

The atmospheric drag on an object will continually decrease the osculating semi-major axisa,according to,22

where μ is the product of the gravitational constant and the mass of the Earth,vis the speed of the object,evis the unit vector in the direction ofv,and˙vDis the acceleration of the object due to drag,given by g

where ρ is the atmospheric density at that altitude,Vis the atmospheric wind velocity vector,andev-Vis the unit vector of the object’s motion relative to the atmospheric wind.The ballistic coefficients BC is defined as

whereCDis the drag coefficients of the object,Ais the cross sectional area of the object in the direction of the object’s motion relative to the atmosphere,andmis the object’s mass.

Based on Eq.(1),the change rate of the mean semi-major axis due to atmospheric drag is derived by Picone et al.17as

whereamis the mean semi-major axis,and the dimensionless wind factorFis given by

The wind factorFhas a good approximation,given by King-Hele22in the form of

wherereis the distance of the object from the center of the Earth,neis the Earth rotation rate andiis the angle of inclination of the orbit.The perigee distancerp0and perigee velocityvp0can be used in place ofrandv,respectively.

Integrating Eq.(4)from timet1tot2,one has

The numerical integration of the above equation can be written as

where Δtis the time step for the numerical integration.

Thus,the BC can be computed by

2.2.Unstable directly estimated results

The BC value can be directly computed by Eq.(9).The time step Δtis 5 min.Due to the limitation in accuracy ofaD,the integral interval ΔT=ti+1-tishould be suf fi ciently long in order to obtain stable results;The time interval is at least

three days,and the start times of two neighboring time intervals have a span of at least 12 h.

Some parameters used in Eq.(9)are approximate values and contain error.The atmospheric density ρ is difficult to

obtain since the dynamics of the upper atmosphere are not completely understood and the atmosphere density cannot be measured by real-time observation.The conventional way to calculate ρ is to use the value derived from the empirical model,which induces large error.Because of the accuracy of approximated parameters,estimated BC values have unstable oscillations,as shown in Section 4.2.1.However,it can be seen that the estimated BC values of two objects with similar orbital elements,(e.g.altitude,eccentricity and inclination)have similar oscillations.This is because the errors from the approximated parameters and empirical models have similar effects on the objects.With the use of division,ratios between different objects offset the effects of error,hence accurate and stable results are obtained.This process significantly mitigates the effects from the approximation and empirical models.

2.3.Comparison between different objects

The ballistic coefficients ratios can be calculated by comparison of the variation in mean semi-major axes between different objects.The comparison approach between BC values of objectsiandjis given by

This comparison is based on the direct estimation method in Eq.(9)and mitigates error.In order to obtain accurate ratios,two objects should have similar orbital elements.The proof is as follows:

Since the atmospheric density varies with the latitude,and variance is significantwith longitude,the comparison approach can be more specific when considering the former.In order to compare two objects in an analogous atmospheric environment,this approach can be used to compare computed BC values at similar latitude.

One orbital period is divided intonTparts with time step Δt=Suppose there arenPorbital periods in time interval ΔT.The second fraction in Eq.(10)can be written as the sum of different parts at different latitudes,i.e.,

Letrdenote the proportional error,which is the proportion of estimated value from approximation and empirical models to the true value.The above equation becomes

Suppose objectiandjhave similar orbital elementsparts at similar latitude have similar proportional errorThe comparison within the parts at similar latitude is given by

The proportional errorrhas been cancelled out.But,the variation ΔaDat different latitudes cannot be obtained precisely.The comparison has to be made across the whole time interval with variation Δt2t1aD.If the satellites are running in near circular orbit,one has

Then,the effects ofrcan be cancelled out in Eq.(12),i.e.,

In each time interval,an estimate of the ballistic coefficient ratio can be obtained through this comparison approach.The final ratio takes the average of the results obtained in all time intervals.When the BC value of one object is known,the BC value of the other can be calculated by

This process mitigates the error when two objects are at similar orbits.The more similar the orbital elements are,the better the results.

2.4.Two-line element sets

The most widely used orbital data are the Two-line Element sets(TLEs),which are available to general public.The TLEs provides a vast storehouse of standard,accessible and easyto-use orbital data for space research.The TLEs data used here are available on the web:http://www.celestrak.com.The TLE consists of six orbital elements,and these orbital elements are similar but not identical to the classical elements.They are mean orbital elements without the periodic terms.23

The TLE data does not include the semi-major axis item,but has the mean mean-motionnmwith units of revolutions per day(rev/day).The mean mean-motion does not include the long-periodic and short-periodic terms,thus the variation in the mean mean-motion of TLEs is mainly due to the drag effect.The TLE data can be used to estimate atmospheric drag and derive the atmospheric density models.24

The mean mean-motionnmis the Kozai mean value,and the relationship between the mean mean-motion and the mean semi-major axis is written as

2.5.Selection of time interval

The TLEs may be released two or three times per day and typically occur at twelve o’clock a.m.or p.m.The time intervals of two compared objects are chosen between similar release times.

Due to limitations in the accuracy of TLEs,every interval ΔTkis at least three days,and the start time of two neighboring time intervals ΔTk-1and ΔTkshould have at least a 12 h span.Selection of time intervals is shown in Fig.1.

The number of orbital periods in the chosen time interval needs to be calculated.The mean mean-motionnmin TLEs is measured in revolutions per day,thus the integer part is the product ofnmand ΔTrounded to the nearest integer.

If the two TLE sets do not have similar release times over a period,the time interval is shorter.Thus,in order to obtain the variationa fi tting process is necessary and introduced in Section 3.

2.6.Calculation of the parameters

At each time step,the factorsak,Fkand ρkcan be considered constant.The semi-major axisakvaries fromak（t1）toak（t2）during time interval ΔTk,andak（t）can be obtained from the fitted curve ofa.The wind factorFkin each time step can be calculated by Eq.(6).

The atmospheric density ρkis estimated by empirical model.In this work,the naval research laboratory mass spectrometer and incoherent scatter(NRLMSISE)model is chosen.The NRLMSISE model has been revised to include total mass densities from satellite accelerometers and orbit determination and designated the NRLMSISE-00 model.25

In calculating atmospheric density,the solar 10.7 cm radio flux(F10.7)is used to account for variations in solar extreme ultraviolet(EUV),and theapindex is a proxy for geomagnetic disturbances.The measurements of these two parameters also have error,which may be carried to estimated ρ.The data ofF10.7andapfrom 2000 to 2016 are available on the web:www.celestrack.com.

3.Fitting

Due to measurement and fitting error of TLEs,the variation ofamis not monotonic decreasing,but presents irregular shaking and even increases over a short period time.These lead to singular results and are also the reason to use a variation ofamwith a sufficiently long time.Thus,the polynomial fitting presented in Bar-Shalom et al.26is used to f i t the change ofam,and the fitting order should be chosen appropriately.The LS technique is used for polynomial fitting.

3.1.Goodness-of-f i t

The residuals in an LS estimation problem are also called the goodness-of-f i t or fitting error,which is given by

Under Gaussian assumption,the sum of squares of fitting error is chi-square distributed withknz-nxdegrees of freedom,wherekis the number of measurements of dimensionnz,andnxis the dimension of the parameter vector.

The noise covariance matrix of TLEs is unknown,but can be estimated together with the regression coefficients.Also,as a result of the many studies done on the accuracy of TLEs,the noise of TLEs can be considered within 102–103m.27,28

3.2.Underfitting and overfitting

If the order of polynomial is too low,i.e.,underfitting,the f i t will be poor.If the order of the polynomial is too high,i.e.,overfitting,some of the estimated parameters are statistically insignificant and thus become ‘noise”.

(1)Test for underfitting:The order of the polynomial f i t is too low when the fitting error is too high,i.e.

wherecis the tail probability of the chi-square density such that the probability of aknz-nxdegrees of freedom chisquare random variable exceeding it is α(5%).

(2)Test for overfitting:The order of the polynomial is too high when some estimates of coefficients are statistically insignificant.With normal noises,the estimate of the i th component of the parameter vector is given by

whereN（xi，Pii（k））denotes a normal with meanxiand variancePii（k）.The parameter estimate significance test is the test between the two hypothesesH0:xi=0 andH1:xi≠0.ThenH1can be accepted only if

The parameter estimate significant test implies the 95%two-sided probability region.The thresholdc′is the value that the probability of a standard Gaussian random variable exceeding it isIf the test shows the parameter is statistically insignificant,lower the dimension parameter.

A segment of semi-major axis variation of TANSUO-1(28220)over 120 days is used as a fitting example.The results are illustrated in Fig.2 with the 2σ confidence region(95%)shown.

For the second-order and third-order polynomial f i t,the fitting error is large.For the fourth-order case,the fitting error is acceptable and all the parameter estimates are significant.Compared with the fourth-order case,the fifth-order f i t has only slightly improved but the improvement is statistically insignificant,which can be seen in the confidence region widening.This is an overfitting case;thus,for this fitting segment,the fourth-order f i t is accepted.

3.3.Piecewise fitting

In order to achieve the best approximation,piecewise fitting is used.But,the fitting result is not continuous at the boundary point between neighboring pieces.To smooth the boundary,each piece is extended forward and backward to make the boundary point meet the interpolation condition.Then,interpolation algorithms can be used at the boundary of each piece to keep the continuity of the boundary point.29

3.4.Outliers elimination

Consider the measurement and fitting errors in TLEs,there may be outliers.The TLEs with obvious abnormal values are not used for estimation.Moreover,there may be continuous outliers over a time period;the TLE data of this time period are not used for estimation and the ratios with abnormal values are eliminated from final results.

4.Numerical experiments

4.1.Test space objects

Ten satellites are chosen to test the proposed method.The parameters of test objects are known,and they have stable area-to-mass ratios with regular shape(spherical or square).Their orbital data are shown in Table 1.

The ballistic coefficients of these satellites are known,or the satellites have regular shape and their ballistic coefficients can be easily computed.

(1)For ERS-2,its reference BC is computed with a mass of 2500 kg,a cross-sectional area at 21 m2in its motion direction,and a CDvalue of 3.1 as used in Bennett et al.30

(2)The shape of Clementine is nearly spherical with diameter of 1 m.A mass of 50 kg and a CDvalue of 2.2 are used to calculate its reference BC value.

(3)Bowman et al.31give the estimation of BC values for the CHAMP and GRACE satellites in which the CDhave values around 3.3–3.4.

(4)The Starshine-3 is a sphere,and its drag coefficients variability is analyzed in Bowman and Moe.32

(5)Both the TANSUO-1 and NAXING-1 satellites have nearly spherical shapes with diameters of 1 m and 0.5 m,respectively.The mass of TANSUO-1 and NAXING-1 are 204 kg and 25 kg,respectively.The CDuses a value of 2.2 for both.

(6)The SWAS is a 300 kg and 2-m-diameter spherical satellite,and its CDis 2.2.

(7)The RESURS DK-1 is a 6 m×3 m cylindrical satellite with mass of 6804 kg,and its cross-sectional area is 15 m2.Similar to CHAMP and GRACE,a value of 3.3 is used for its CD.

4.2.Results and analyses

The ten test satellites are separated into 5 groups,and each two satellites with similar orbital elements are processed as a satellite pair.Some satellite pairs contain two satellites launched by the same rocket that have either similar BC values(27391 and 27392)or different BC values(28220 and 28221).The other paired satellites have no relationship with each other.Also,two satellites(26405 and 26929)with different inclination are compared to test whether the proposed method can be used for such cases.

4.2.1.Directly estimated ballistic coefficients values

The directly estimated ballistic coefficients of ten objects are shown at the top of Figs.3–7.The directly estimated BC values of two paired satellites are shown in the same figures.Due to error from approximation and empirical models,the directly estimated BC values are unstable and fluctuant in the data period.True values cannot be determined,and the averaged values over a long period of time are taken as estimated BC values and are presented in Table 2.The directly estimated BC values of 27391 and 27392 have an apparent increase in Fig.4,and the mean estimated BC values of these two after 2010 are 0.00626 and 0.00627,respectively.These values aremore accurate and also indicate that the directly estimated results are unstable.

Table 1 Data of test LEO objects.

4.2.2.Ballistic coefficients values computed by ratios between different objects

The ratios of BC values in each satellite pair and the corresponding probability density functions(PDF)are shown in the bottom left and bottom right of Figs.3–7,respectively.From these figures,it can be seen that the directly estimated ballistic coefficients show similar fluctuations to those space objects with similar orbital elements;the ratios alleviate the fluctuations and stable results are obtained.The variation in relation of parameters between different objects with similar orbital elements is also confirmed.The estimated BC values by ratios of the ten objects are presented in Table 3.Compared with the directly estimated results in Table 2,those in Table 3 demonstrate that the proposed method yields accurate estimates and outperforms the direct estimation method.

Table 2 Directly estimated BC values of ten test objects.

5.Experiments for LEO objects with higher altitude

5.1.Test LEO objects higher than 600 km in perigee altitude

The proposed method has been validated for LEO objects lower than 600 km in perigee altitude.Some LEO objects with perigee altitude higher than 600 km are also used to test the proposed method;this data are presented in Table 4.

Two geodetic satellites,Stella and Westpac,are chosen as test objects.Stella and Westpac are both spherical objects with a diameter of 240 mm with masses of 48 and 23 kg,respectively.Harrison and Swinerd33analyzed the drag coefficients of Stella by using gas-surface interaction assumptions.Bennett et al.30computed the BC value of Westpac withCD=2.3.

Table 3 Estimated BC values using ratios between different objects.

Table 4 Data of test LEO objects higher than 600 km.

Two amateur satellites are also chosen due to their regular shape and easy to calculate BC values.The main bodies of Eyesat1 and Itamsat are cubic objects with side lengths of 230 mm,and the masses of these two objects are 11.8 and 11.2 kg,respectively.TheCDvalue of 2.3 and cross-sectional area at 0.1 m2are used for both.For more details on the amateur satellites refer to:www.amsat.org.

5.2.Proportional errors are cancelled out in comparisons

For space objects with perigee heights of about 800 km,solar radiation pressure may be at the same level as atmospheric drag force.But,over long timescales(e.g.,annual),this pressure can integrate to a small number.This can induce variations in semi-major axes on timescales of 30 days to one year.17Thus,for short-term comparison,the effect of solar radiation pressure is small.

The acceleration of the solar-radiation force is given by whereCRis the reflectsivity,pSR=4.51×10-6N/m2is the solar-radiation pressure,A′is the exposed area to the Sun,ande⊙is the unit vector of the satellite relative to the Sun.

This equation has an itemCRsimilar to the ballistic coef-that is to say the acceleration of solarradiation force is proportional to the area-to-mass ratioWhile the space object has a regular shape,e.g.sphere,it can be assumed that the frontal areaAis the same as the area exposed to the SunA′.Thus,the acceleration of drag force and the solar-radiation force can be considered proportional to the area-to-mass ratio,and the variation of semi-major axiscan still be considered proportional to the area-tomass ratioThus,the results indicate the proposed method also can be used for LEO objects higher than 600 km in perigee altitude.

Table 5 Directly estimated BC values of test LEO objects higher than 600 km.

LetCR=sCD,and assumeA=A′.Regardless of the direction of the acceleration,the acceleration of the atmospheric drag force is given by

Thus,the ratio of the acceleration of two objects can be written as

When the comparison time span is short,the parameters in Eq.(24)can be considered constant.When the orbits of two objects are close enough,the first ratio in Eq.(24)is approximate to 1,and the effects of solar-radiation pressure are cancelled out.Thus,when two objects are in a similar perigee altitude higher than 600 km,Eq.(10)is still valid and can be used to compute BC ratios.

5.3.Results and analyses

Numerical experiments are used to validate the proposed method for objects higher than 600 km in perigee altitude.The directly estimated BC values for the four objects are presented in Table 5.Due to influence from other forces,computed BC values should consider more than only the atmospheric drag force.The estimated BC values presented in Table 5 show large relative errors.

In order to test the proposed method,each of the four satellites takes turns being the reference object to estimate the BC values of other three.The estimated BC values from using these four satellites are shown in Tables 6–9.It can be seen that values estimated by ratios have high accuracy,which indicates good performance of the proposed method for objects with perigee altitude higher than 600 km.

6.Experiments for space debris

Debris objects originating from one source can be defined together as a debris group.Since the debris in a group has similar orbital elements,stable BC ratios can be obtained between them.Three pieces of debris are used as research objects.Their data are presented in Table 10.Their directly estimated BC values and the PDF of the BC ratios are shown at the top and bottom of Fig.8,respectively.

Table 6 Estimated BC values using 22824.

Table 7 Estimated BC values using 22825.

Table 8 Estimated BC values using 22826.

Table 9 Estimated BC values using 25398.

Table 10 Data of test debris objects.

From Fig.8,it can be seen that the ratios of BC values between debris objects in this debris group are stable,due to the similar orbital elements and space environment of debris objects.This result confirms that these three objects originated from the same source and can be grouped together.This process gives information for the debris group that is a new resource and gives theoretical and methodological support for the catalogue of debris group.Based on this information,unknown debris can be determined to belong to a particular debris group.

7.Discussion

In order to verify estimation results,the ballistic coefficients ratios can be used inversely for prediction.Also,this process can be used to predict the secular change in semi-major axes of space debris without propagators when they are not observed for a period of time.In time intervalt1tot2,if the ballistic coefficient ratioRand the variation of the mean semi-major axisof objectjare known,the variation of the semi-major axisof objectican be predicted by

In order to demonstrate the results clearly,the probability density of ballistic coefficients ratios is used to show outcomes with different probabilities.The Sequential Monte Carlo Implementation is used to propagate the weighted particles that approximate the probability density.34In this process,the probability density can be constructed using particle sampling.A number of particles are sampled,uniform from the probability density of the ballistic coefficients ratio,and the probabilities are indicated by the weights of particles.The particles are propagated through prediction equation,Eq.(25).

The satellite pairs used to compute the ballistic coefficients ratios are used here for prediction.One satellite pair is used as a reference object to predict the semi-major axis of another.The satellites 25978,27392,28221 and 25560 are used as the reference objects to predict the variation in the semi-major axes of satellites 23560,27391,28220 and 29228,respectively.The two-year prediction results from 2014 to 2016 are shown in Fig.9.

Since the ballistic coefficients ratio is propagated as probability density,the deeper the gray,the higher the probability of the space debris appearing at that place.The red line is the true change ofa.It can be seen that the long-term change ofapredicted by inversely using the proposed approach has good results without an orbital propagator model.

8.Conclusions

Due to the rapidly growing amount of space debris,its research is of great importance.Ballistic coefficients are a sig-nificant parameter for the cataloguing,prediction and mitigation of space debris.A new method is proposed in this paper to estimate the ratios of ballistic coefficients between different space objects.The BC values of research objects can be estimated when the parameters of reference objects are known.In addition,ballistic coefficients ratios are used in the research of debris grouping.The ratios provide new information and give theoretical and methodological support for debris group cataloguing.This method has been verified through tested of objects with known parameters.The results validate the proposed method for LEO objects both lower and higher than 600 km in perigee altitude.Compared with direct estimation methods,estimation by ratios demonstrates advantages in stability and accuracy since it mitigates error from approximating parameters of both the compared objects.In order to verify estimated BC values,this approach is used inversely to predict secular change of semi-major axes.The estimation by ratios approach provides a new way to predict long-term change of semi-major axes without an orbital propagator model.

Acknowledgements

The authors would like to gratefully acknowledge the research support from Applied Astronomy Research Group,Yunnan Observatories,Chinese Academy of Sciences and the grant support from the National Natural Science Foundation of China(No.61372162).

1.Kessler DJ,Cour-Palais BG.Collision frequency of artificial satellites:the creation of a debris belt.J Geophys Res:Space Phys1978;83(A6):2637–46.

2.Su SY,Kessler D.Contribution of explosion and future collision fragments to the orbital debris environment.Adv Space Res1985;5(2):25–34.

3.Kessler DJ.Orbital debris environment for spacecraft in low earth orbit.J Spacecraft Rock1991;28(3):347–51.

4.Klinkrad H,Jehn R.The space-debris environment of the earth.ESA J1992;16:1–11.

5.Rossi A,Anselmo L,Cordelli A,Farinella P,Pardini C.Modelling the evolution of the space debris population.Planet Space Sci1998;46(11):1583–96.

6.Petro AJ.Techniques for orbital debris control.J Spacecraft Rock1992;29(2):260–3.

7.Wu Z,Hu R,Qu X,Wang X,Wu Z.Space debris reentry analysis methods and tools.Chin J Aeronaut2011;24(4):387–95.

8.Shen S,Jin X,Hao C.Cleaning space debris with a space-based laser system.Chin J Aeronaut2014;27(4):805–11.

9.Ishige Y,Kawamoto S,Kibe S.Study on electrodynamic tether system for space debris removal.ActaAstronaut2004;55(11):917–29.

10.Huang J,Hu W.MCMC-particle-based group tracking of space objects within Bayesian framework.Adv Space Res2014;53(2):280–94.

11.Hoots FR,Roehrich RL.Space track report no.3 models for propagation of NORAD element sets.Peterson:Aerospace Defense Command,United States Air Force.1980.p.1–79.

12.Vallado DA,Crawford P,Hujsak R,Kelso TS.Revisiting spacetrack report#3AIAA/AAS astrodynamics specialist conferenceandexhibit,2006August21–24;Keystone,Colorado,USA.Reston:AIAA;2006.p.1–88.

13.Schumacher Jr P,Glover RA.Analytic orbit model for US naval space surveillance:an overview.Proceedings of U.S.-Russian second space surveillance workshop;1996 July 4–6;Poznan.1996.p.684 66–105.

14.Linares R,Jah MK,Crassidis JL.Space object area-to-mass ratio estimation using multiple model approaches.Adv Astronaut Sci2012;144:55–72.

15.Bowman BR.True satellite ballistic coefficients determination for HASDM.AIAA/AAS astrodynamics specialist conference;2002 August 5–8;Monterey.Reston:AIAA;2002.p.774–93.

16.Storz MF,Bowman BR,Branson MJI,Casali SJ,Tobiska WK.High accuracy satellite drag model(HASDM).Adv Space Res2005;36(12):2497–505.

17.Picone J,Emmert J,Lean J.Thermospheric densities derived from spacecraft orbits:accurate processing of two-line element sets.J Geophys Res:Space Phys2005;110(A13):103–15.

18.Saunders A,Swinerd GG,Lewis HG.Deriving accurate satellite ballistic coefficients from two-line element data.J Spacecraft Rock2012;49(1):175–84.

19.Sang J,Bennett JC,Smith CH.Estimation of ballistic coefficients of low altitude debris objects from historical two line elements.Adv Space Res2013;52(1):117–24.

20.Montenbruck O,Gill E.Satellite orbits:models,methods and applications.Germany:Springer Verlag;2012.p.83–5.

21.Liu L.Orbit theory of spacecraft.Beijing:National Defense Industry Press;2000,p.114–9[Chinese].

22.King-Hele D.Satellite orbits in an atmosphere:theory and applications.Glasgow:Blackie;1987.p.155–8.

23.Vallado DA.Fundamentals of astrodynamics and applications.4th ed.Hawthorne:Microcosm;2013.p.140–2.

24.Emmert J,Picone J,Meier R.Thermospheric global average density trends,1967–2007,derived from orbits of 5000 near-Earth objects.Geophys Res Lett2008;35(5):84–102.

25.Picone J,Hedin A,Drob DP,Aikin AC.NRLMSISE-00 empirical model of the atmosphere:statistical comparisons and scientific issues.J Geophys Res:Space Phys2002;107(A12):15–6.

26.Bar-Shalom Y,Li XR,Kirubarajan T.Estimation with applications to tracking and navigation.New York:Wiley;2001,p.145–61.

27.Flohrer T,Krag H,Klinkrad H.Assessment and categorization of TLE orbit errors for the US SSN catalogue.Advanced maui optical and space surveillance technologies conference.2008.p.513–24.

28.Levit C,Marshall W.Improved orbit predictions using two-line elements.Adv Space Res2011;47(7):107–15.

29.Shi HL,Yan YH.Extended interpolation method and its applications in piecewise approximations.Comput Appl Math1992;12:229–36.

30.Bennett JC,Sang J,Smith CH,Zhang K.Accurate orbit predictions for debris orbit manoeuvre using ground-based lasers.Adv Space Res2013;52(11):1876–87.

31.Bowman BR,Marcos FA,Moe K,Moe MM.Determination of drag coefficients values for CHAMP and GRACE satellites using orbit drag analysis.Adv Astronaut Sci2008;129(1):147–66.

32.Bowman BR,Moe K.Drag coefficients variability at 175–500 km from the orbit decay analyses of spheres.Adv Astronaut Sci2008;129(1):207–22.

33.Harrison I,Swinerd G.A free molecule aerodynamic investigation using multiple satellite analysis.Planet Space Sci1996;44(2):171–80.

34.Doucet A,de Freitas N,Gordon N.Sequential Monte Carlo methods in practice.New York:Springer Verlag;2001.p.17–33.

5 September 2016;revised 28 October 2016;accepted 27 November 2016