Abstract
The paper aims at discussing techniques for administering one implementation issue that often arises in the application of particle filters: sample impoverishment. Dealing with such problem can significantly improve the performance of particle filters and can make the difference between success and failure. Sample impoverishment occurs because of the reduction in the number of truly distinct sample values. A simple solution can be to increase the number of particles, which can quickly lead to unreasonable computational demands, which only delays the inevitable sample impoverishment. There are more intelligent ways of dealing with this problem, such as roughening and prior editing, procedures to be discussed herein. The nonlinear particle filter is based on the bootstrap filter for implementing recursive Bayesian filters. The application consists of determining the orbit of an artificial satellite using real data from the GPS receivers. The standard differential equations describing the orbital motion and the GPS measurements equations are adapted for the nonlinear particle filter, so that the bootstrap algorithm is also used for estimating the orbital state. The evaluation will be done through convergence speed and computational implementation complexity, comparing the bootstrap algorithm results obtained for each technique that deals with sample impoverishment.
1. Introduction
The orbit of an artificial satellite is determined using real data from the Global Positioning System (GPS) receivers. In the orbit determination process of artificial satellites, the nature of the dynamic system and the measurements equations are nonlinear. As a result, it is necessary to manage a fully nonlinear problem in which the disturbing forces as well as the measurements are not easily modelled. In this orbit determination problem, the variables that completely specify a satellite trajectory in the space are estimated, with the processing of a set of pseudorange measurements related to the body.
A spaceborne GPS receiver is a powerful resource to determine orbits of artificial Earth satellites by providing many redundant measurements, which ultimately yields high degree of the observability to the problem. The Jason satellite is a nice example of using GPS for space positioning. Through an on-board GPS receiver, the pseudoranges (error corrupted distance from satellite to each of the tracked GPS satellites) can be measured and used to estimate the full orbital state.
The bootstrap filter is a particle filter whose central idea is to express the required probability density function (PDF) as a set of random samples, instead of a function over state space [1–3].
Numerous strategies have been developed for solving the particles degeneracy (or sample impoverishment) problem that often arises in particle filter applications like introduction of a risk-sensitive particle filter as an alternative approach to mitigate sample impoverishment based on constructing explicit risk functions from a general class of factorizable functions [4]; incorporation of genetic algorithms into a particle filter [5, 6]; and many others [7–9]. All these strategies, although extremely interesting and suitable for the orbit determination problem, are not in the scope of this work. Here, the option was done for studying two classical methods to solve (or try to solve) the degeneracy problem: roughening and prior editing.
Herein, the main goal is to analyze the bootstrap filter behavior for the highly nonlinear orbit determination problem. Its simulation results are compared taking into account the sample impoverishment. A reference solution is a bootstrap particle filter (BPF) applied to orbit determination that has already been compared to the unscented Kalman filter solution for the same problem and works well for the analysis of the sample impoverishment issue [10].
2. Particle Filter
The particle filter was designed to numerically implement the Bayesian estimator [2]. The Bayesian approach consists of constructing the PDF of the state based on all the available information, and, for nonlinear or non-Gaussian problem, the required PDF has no closed form. The bootstrap filter represents the required PDF as a set of random samples, which works as an alternative to the function over state space. This filter is a recursive algorithm for propagation and update of these samples for the discrete time problem. The Bayes rule, the key update stage of the method, is implemented as a weighted bootstrap [1].
The main idea of the BPF is intuitive and direct. At the beginning, particles () are randomly generated, based on the known initial PDF . At each step of time , the particles are propagated to the next step using the dynamics equation [2]. After receiving the measurement at time , the PDF is evaluated. That is, the conditional relative likelihood of each particle is calculated. If an -dimensional measurement equation is given as and is a Gaussian random variable with a mean of zero and a variance of , , then a relative likelihood that the measurement is equal to a specific measurement , given the premise that is equal to the particle , can be computed as follows [2]:
In (1), the symbol means that the probability is directly proportional to the right side. So if the equation is used for all particles , then the relative likelihood that the state is equal to each particle is correct. The relative likelihood values are normalized to ensure that the sum of all likelihood values is equal to one. Next, a new set of randomly generated particles is computed from the relative likelihood . In the resampling step, roughening was used, in order to prevent sample impoverishment. At this point, there is a set of particles that are distributed according to the PDF , and any desired statistical measure of it can be computed [2].
The particle filter, adjusted to the orbit determination problem, can be summarized as follows. (1) The dynamic and the measurement equations are given as where and are independent white noise processes with known PDFs. (2) initial particles () are randomly generated on the basis of the known initial state PDF . is a parameter chosen as a trade-off between computational cost and estimation accuracy [2]. (3) For , (a) in the time propagation step, obtain a priori (predicted) particles , using the dynamics equation and the PDF of the process noise, both known: where each noise vector, , is randomly generated on the basis of the known PDF of ; (b) compute the relative likelihood of each particle , conditioned on the measurement , using the nonlinear measurement equation and the PDF of the measurement noise, as in (1); (c) normalize the relative likelihood values: (d) in the resampling step, generate a set of a posteriori (resampled) particles , on the basis of the relative likelihood ; (e) now, there is a set of particles distributed according to the PDF , and mean and covariance statistical measures can be computed.
In the implementation of the bootstrap filter, there is only a small overlap between the prior and the likelihood.
There are some procedures that may be implemented for combating the consequent reduction in the number of truly distinct sample values, such as increasing the number of particles, roughening, and prior editing [1]. Here, they were implemented: a bootstrap particle filter with resampling (PF); a PF with roughening (PFR); and a PFR with prior editing (PFPE), in order to evaluate roughening and prior editing strategies for dealing with sample impoverishment.
2.1. Roughening
Roughening will be the first remedy for sample impoverishment to be discussed. It restrains the resampled particles spread (a posteriori particles) by adding random noise to them, which is similar to adding artificial process noise to the Kalman filter [2]. In roughening approach, the a posteriori particles are modified, after the resampling step, as follows: is a zero-mean random variable (usually Gaussian); is a constant tuning parameter; is the number of particles; is the state space dimension; and is a vector of the maximum difference between the particle elements before roughening. The th element of the vector is given as where is the step time and and are particle numbers.
The tuning parameter choice is a compromise. Being too large, a value would blur the distribution, but being too small, it would produce tight clusters of points around the original particles [1]. In this paper, .
2.2. Prior Editing
Prior editing can be tried if roughening does not prevent sample impoverishment. Such approach edits the a posteriori particles from the prior time instant, (after roughening), if the a priori particle from actual instant, , does not satisfy a coarse, pragmatic acceptance test [1]. Therefore, this procedure artificially boosts the number of samples of the prior editing in the neighborhood of the likelihood, for if an a priori particle is in a region of state space with small , it is rejected. Then, the a priori rejected particle can be roughened as many times as required, according to (5), until it is in a region of significant [2]. The prior editing was implemented as follows [1]:(a)Pass the resampled sample from previous instant, , through roughening and system model to generate the predicted sample from current instant, .(b)Calculate , the residual between the true and the predicted measurements, for the th particle of the sample, considering that the actual instant observation is available.(c)If the magnitude of is higher than six standard deviations of the measurement noise, then it is highly unlikely that is chosen as an a posteriori particle. In this case, is rejected, and is roughened again and passes one more time through dynamic model to generate a new a priori particle . As has already passed through roughening and generated a rejected predicted particle, this procedure may be repeated while is in a region of no negligible probability.
Due to the high computational cost involving prior editing, such approach was done only once. It is important to make it clear that, here, the th particle is, in fact, an -dimensional vector, while a sample is a matrix, where the th state variable is represented by particles.
The accommodation of roughening and prior editing in the bootstrap particle filter algorithm can be schematized as Figure 1 shows.
3. Orbit Determination
The orbit determination is a process for obtaining values of the parameters that completely specify the motion of an orbiting body (as an artificial satellite), based on a set of observations of the body. It involves nonlinear dynamical and nonlinear measurement systems, which depends on the tracking system and the estimation technique [11, 12]. The dynamical system model consists of describing satellite orbital motion, which includes Earth’s rotation effects and perturbation models and measurements models. These models depend on the state variables initial conditions, as well as a variety of parameters which affect both the dynamic motions as the measurement process [13]. Due to the complexity of the applied models, usually it is not possible to solve such models equations directly for any of these parameters from a given set of observations.
The observation may be obtained from the ground station networks using laser, radar, Doppler, or space navigation systems, as the GPS. The choice of the tracking system depends on a compromise between the goals of the mission and the available tools. In the case of the GPS, the advantages are global coverage, high precision, low cost, and autonomous navigation resources. The GPS may provide orbit determination with accuracy at least as good as the methods using ground tracking networks. The latter provides standard precision around tens of meters and the former can provide precision as tight as some centimetres. The GPS provides, at a given instant, a set of many redundant measurements, which makes the orbit position observable geometrically.
After some advances of technology, the single frequency GPS receivers provide a good basis to achieve fair precision at relatively low cost, still attaining the accuracy requirements of the mission operation. The GPS allows the receiver to determine its position and time, geometrically, anywhere at any instant, with data from at least four satellites. The principle of navigation by satellites is based on sending signals and data from the GPS satellites to a receiver located on board the satellite whose orbit needs to be determined. This receiver measures the travel time of the signal and then calculates the distance between the receiver and the GPS satellite. Those measurements of distances are called pseudoranges.
The instantaneous orbit determination using GPS satellites is based on the geometric method. In such method, the observer knows the set of GPS satellites position in a reference frame, obtaining its own position in the same reference frame.
3.1. Dynamics Model
In the case of orbit determination via GPS, the ordinary differential equations which represent the dynamic model are, in its simplest form, given traditionally as follows:wherein the variables are placed in the inertial reference frame. In (7), is the vector of the position components (); is the velocity vector; represents the modelled perturbing accelerations; is the white noise vector with covariance ; is the user satellite GPS clock bias; is the user satellite GPS clock drift; and is the noise associated with the GPS clock. The GPS receiver clock offset was not taken into account, so as not to obscure the conclusions drawn in this paper due to introduction of clock offset models in the filters. Indeed, the receiver clock offset was beforehand obtained and used to correct the GPS measurements, so that the measurements are free from the error derived from receiver clock offset.
3.2. Forces Model
There are gravitational and nongravitational forces that affect the orbit of an Earth’s artificial satellite. The main disturbing forces of gravitational nature are the nonuniform distribution of Earth’s mass; ocean and terrestrial tides; and the gravitational attraction of the Sun and the Moon. And the principal nongravitational effects are Earth atmospheric drag; direct and reflected solar radiation pressure; electric drag; emissivity effects; relativistic effects; and meteorites impacts.
The disturbing effects are, in general, included according to the physical situation presented and to the accuracy that is intended for the orbit determination. Here, only a minimum set of perturbations was included which enable analyzing the performance of the particle filter [14]: geopotential [15]; direct solar radiation pressure [16, 17]; and third body point mass effect of the Sun and the Moon [18, 19].
3.3. Observations Model
The nonlinear equation of the observation model iswhere, at instant , is the vector of observations; is the nonlinear function of state , with dimension ; and is the observation errors vector, with dimension and covariance . For the present application, the ion-free pseudorange measurements from the Jason-2 GPS receiver are used. Also, the receiver clock offset was computed before and used to correct the pseudorange measurements. Additionally, the nonlinear pseudorange measurement was modelled according to [20].
4. Results
The tests and the analysis for the bootstrap particle filter and two procedures for avoiding sample impoverishment (roughening and prior editing) are presented. To validate and to analyze the methods, real GPS data from the Jason-2 satellite are used. Ocean Surface Topography Mission (OSTM)/Jason-2 is a follow-on altimetry mission to the very successful TOPEX/Poseidon mission and Jason-1. It is a joint mission between NASA and CNES (French space agency), launched June 20, 2008. Jason-2 has a repeat period of approximately 10 days with 254 passes per cycle. Its nodal period is 6,745.72 sec (near 1.87 hours). Sometimes there may be anomalous or missing data. Occasionally Jason-2 must perform maneuvers to maintain orbit. When the satellite detects something abnormal, it will go into safe hold and will turn off all instruments and no data will be collected [21].
The filters estimated position and velocity are compared with Jason-2 precise orbit ephemeris (POE) from JPL/NASA. The test conditions consider real ion-free pseudorange data, collected by the GPS receiver on-board Jason-2, on October 22, 2010, presenting up to 12 GPS satellites tracked. The tests were limited to 5.5 hours of GPS data spaced 10 s. After that, there was an undesirable data gap which could spoil the test case. Such 5.5-hour arc (near 3 Jason-2 orbital periods) was considered sufficient for evaluating the bootstrap particle filter and roughening and prior editing approaches, in this orbit determination application.
The force model comprises perturbations due to geopotential up to order and degree ; direct solar radiation pressure; and Sun-Moon gravitational attraction. This model of forces is suitable for implementation in orbit determination because of the low computational cost added compared to the improvement in the results accuracy [22–24]. The pseudorange measurements were corrected to the first order with respect to ionosphere.
This work is not a search for results accuracy. It aims at analyzing the application of a bootstrap filter to the orbit determination problem. The analysis is based on comparing the errors in position (translated to the orbital radial, normal, and along-track RNT components) among three solutions:(1)The bootstrap particle filter with resampling, applied without any sample impoverishment procedure (named “PF” in the results).(2)The bootstrap particle filter applied with roughening (named “PFR” in the results).(3)The bootstrap particle filter applied with roughening and prior editing (named “PFPE”in the results).
The RNT system interpretation is straightforward: the radial component “R” points to the nadir direction; the normal “N” is perpendicular to orbital plane; and the transversal (along-track) “T” is orthogonal to “R” and “N,” and it is also the velocity component. Thus, it is possible to analyze what happens with the orbital RNT components and the orbit evolution as well. There is also interest for processing time, in order to analyze the computational efforts face to the accuracy achieved by each algorithm.
Regarding time of processing, , it was expected that the time required for the bootstrap particle filter algorithm was increasing in two scenarios: if the number of particles rises and if roughening and prior editing were added to the algorithm. According to the estimator nature, each element of the state was replaced by an array of particles, where is a trade-off between computational cost and estimation accuracy. Here, tests were done for 17, 100, 300, and 700 particles.
As said before, prior editing recomputes any particle that does not match a specific criterion. Therefore, a relevant test is to observe the instant when each algorithm starts rejecting particles which will be presented. The goal is to verify whether, as a procedure is included, it delays the rejection process. It will also be analyzed whether the number of particles affects the particle rejection, that is, whether its increase may work as an approach for avoiding sample impoverishment.
If is set (so the analysis is per line in Table 1), it is noticeable that , CPU time, measured for PF and PFR is very similar, with a maximum 3% of difference. However, the algorithm that includes prior editing, PFPE, is significantly more costly, reaching 39% of increase. This was expected for PFPE, because of the particles that do not pass the acceptance test and need to return to prior instant. Regarding the effects of increasing the number of particle, for the same algorithm (so the analysis is per column in Table 1), the raise in processing time was 83% from 17 to 100 particles, 67% from 100 to 300, and 59% from 300 to 700. Considering that, for 17, 100, 300, and 700 particles, the increase in is 6, 3, and 2.3 times, respectively; then, increase is not directly proportional to the number of particles used. Therefore, the CPU time and the increase in the number of particles show that the time of processing is related to the chosen number of particles only when the number of particles changes, but it has no relation with the estimator implemented.
For computing time of processing and for all the simulations shown, a computer Intel Core i5-2430M processor of 2.40 GHz, with 2.70 GB of RAM, was used. The program was coded in FORTRAN 77 within operating system Windows XP and compiler Compaq Visual Fortran version 6.1.
When some particles do not reach prior editing criterion (i.e., the magnitude of is higher than six standard deviation of the measurement noise), they need to be edited in the prior instant of time. The first instant when any particle needed to be edited, , was computed, for each algorithm and , as can be verified in Table 2.
In Table 2, the rejection of particles was first detected in the PF algorithm, despite the number of particles used. This was expected, since no procedure in order to combat the number of truly distinct samples reduction was used in this algorithm. And the instant when the first rejection occurred in PFR and PFPE concurred for all changing in the number of particles. This is considerably consistent, since prior editing procedure depends first on roughening implementation. Regarding the increase in the number of particles (from 17 to 100), it is noticeable that, when procedures for avoiding sample impoverishment are adopted (in PFR and PFPE cases), the instant when the first rejection of particles is detected is delayed in 42%. Nevertheless, for the other cases analyzed (300 and 700 particles), such instant remains the same. This suggests that increasing as an approach for minimizing the sample impoverishment issue is a little efficient, with a very high computational burden, as seen in Table 1. For PF, the results were not conclusive. It seems that the higher the number of particles, the faster their impoverishment if nothing is done for avoiding sample impoverishment.
As said before, the number of particles is chosen as a trade-off between computational cost and estimation accuracy. The results in Table 3 present mean and standard deviation of the errors in RNT components evaluated for ; 100; 300; and 700. If only these statistics are analyzed, it is clear that the estimation accuracy improves as the number of particles increases. The largest standard deviation occurred for PFPE (). In the along-track component, a divergence occurrence in all algorithms () was detected, which can be verified by the high mean and standard deviation values obtained in the three cases. Such divergence disappears when is increased. When , the statistics and the errors behavior were very similar in the results obtained by PFR and PFPE versions. By the information presented, it is possible to conclude that if the number of particles is very small, any sample impoverishment avoidance procedure will not be effective. Even more, the PFPE approach as used here (particles edition executed only once) does not present any significant improvement facing the PFR procedure, no matter the number of particles used. And taking into account the higher computational cost for only one edition of particles, if the PFPE is implemented as many times as necessary, the computation burden is enough to contraindicate this procedure use, even if the results are improved.
In order to show the behavior of the errors in terms of RNT coordinates, Figures 2, 3, and 4 are presented. Figure 2 shows the worst result, the solution obtained for PFPE algorithm considering , where the divergence in along-track coordinate is shown. Figure 3 shows a PFR solution for , while Figure 4 presents the solution. It was chosen to introduce PFR solutions because this is the algorithm with better performance during orbit determination process. The errors decrease considerably in all implementations, for 100 or more particles (see Table 3), despite an incongruous behavior close to 1 and 5 hours of processing, from which the filters recover later. This anomalous behavior is detected in all simulations results. In the graphics, blue curves correspond to radial component (R); red to normal (N); and green to along-track (T). Increasing from 300 to 700 was more meaningful to along-track component, as its mean and standard deviation had higher improvement (decreasing in values) and less significant to normal coordinate.
According to Figures 3 and 4, the results for are as competitive as for . And if computational burden of implementing 700 particles is taken into account, 300 particles are sufficient for the evaluation proposed in this paper.
Despite the undesirable data gap near 5.5 hours of GPS data, in order to properly evaluate the fittest algorithm (PFR) regarding number of particles, two other graphics were generated. Figures 5 and 6 present and , the errors in position and velocity, respectively, for 24 hours of implementation. The results obtained for (green curves) are compared with (dark blue curves). According to the results, it is clear that a higher number of particles are important for improving the results, since the amplitude of errors diminished considerably when was used. And despite the troublesome data, for 300 particles case, PFR can recover and continue the estimation process, without divergence, which is an indication of the particle filter robustness. In each figure, and represent, respectively, the absolute value of the errors in position and in velocity, in the ECEF (Earth-Centered Earth-Fixed) reference frame, also known as ECR (Earth-Centered Rotational).
5. Conclusions
A Bayesian bootstrap particle filter was applied for the satellite orbit determination problem, using a set of GPS measurements. The development was evaluated taking into account performance and computational burden. The bootstrap filter results, implemented with resampling, were compared with two other versions, which include roughening and prior editing, aiming at avoiding sample impoverishment.
With regard to time of processing, results showed that PFPE algorithm requires greater time than PF, which is as competitive as PFR. Since the number of particles is chosen as a trade-off between computational cost and estimation accuracy, the CPU time is related to the variation in number of particles chosen.
When the number of particles used is analyzed, it is also possible to conclude that if it is very small, any sample impoverishment avoidance procedure will not be effectual. And since PFPE approach executes particles edition only once here, it does not present any significant improvement facing the PFR procedure, no matter the number of particles used. Additionally, if the computational costs between PFR and PFPE (as implemented) are compared, the computation burden is enough to contraindicate PFPE implementation as many times as necessary, even if the results are improved.
Results confirm that the greater the number of particles, the better the estimation accuracy. The best result was achieved for , in the three versions of the estimator, although a higher computational effort was demanded in this case. Therefore, when results are compared, it is possible to assure that 300 particles are enough to achieve the accuracy level aimed in this paper.
In order to obtain a better bootstrap particle filter performance, especially in terms of estimation accuracy, adjustments in the many filter variants might be done for improving its efficiency. Such adjustments are directly related to the knowledge about the filter. Other strategies can also be tried for solving implementation issues such as sample impoverishment. Another approach for improving particle filtering is to combine it with another filter such as extended Kalman filter or unscented Kalman filter. In this approach, each particle is updated at the measurement time using the extended or the unscented filter, and then resampling is performed using the measurements.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors wish to express their consideration to the National Institute for Space Research (INPE) that kindly provided everything necessary for this paper to be developed. The authors are also grateful to São Paulo Research Foundation (FAPESP), for the support under Contract no. 2012/21023-6, and to National Council for Scientific and Technological Development (CNPq), for the support under Contract no. 303119/2010-1.