ANALYZING METEOROID FLIGHTS USING PARTICLE FILTERS

A meteoroid is a small object moving in interplanetary space. When one enters Earth's atmosphere, it creates a bright phenomenon called a meteor, fireball,  or bolide (depending on brightness). The interaction of this material with our atmosphere provides us with an opportunity to observe and study a portion of interplanetary material that would otherwise be inaccessible to us. Telescopes cannot image millimeter- to meter-sized objects, and discoveries of asteroids tens of meters in size constitute a tiny fraction of the predicted population (Harris 2012). Determining the physical state of this material in our atmosphere, its strength and mass distribution, and its velocity frequency distribution provides a unique window on cometary and asteroidal material in the inner solar system. In order to derive that data, we need to model the meteoroid–atmosphere interaction.

A set of idealized equations govern how a single meteoroid body will respond in terms of velocity and mass loss. The amount of deceleration experienced by a meteoroid is related to its shape and bulk density via a shape density parameter, $\kappa ={c}_{d}A/2{\rho }_{m}^{2/3}$, where c_d is the aerodynamic drag coefficient,³ A is the shape parameter as described by Bronshten (1983), and ρ_m is the bulk density of the meteoroid. Both ablation and gross fragmentation of the meteoroid are responsible for loss of mass. Gross fragmentation is hard to predict and is linked to the strength of the object. Ablation can be quantified through the ablation parameter—σ, which is defined as $\sigma =\tfrac{{c}_{h}}{{H}^{* }{c}_{d}}$ (see footnote 1) (where c_h is the coefficient of heat and H* the enthalpy of vaporization).

If the meteoroid survives this luminous trajectory or bright flight, there is the possibility of recovering a meteorite on the ground. Dedicated fireball camera networks such as the Desert Fireball Network (DFN) in Australia (Bland et al. 2012) allow triangulated trajectories of larger meteoroid bodies to be observed. Special shutters are used (in the case of the DFN, a liquid crystal shutter using modulated sequences; Howie et al. 2016) to encode timing throughout the trajectory. Being able to predict the final state of the meteoroid is paramount to determining whether there is any recoverable material, and is a necessary input to so-called dark flight modeling (the process by which data from the luminous trajectory are converted into a fall line on the ground using atmospheric wind models), enabling likely search areas to be defined (Ceplecha 1987). Accurately calculating a trajectory also allows the orbit for that body to be determined. Meteorites with orbits are rare: less than 0.05% of all meteorites. Knowing a meteorite's pre-atmospheric orbit gives contextual information to the picture they provide on early solar system formation. Over time, the statistical analysis of calculated orbits may also assist in planetary defense of asteroidal debris streams.

Determining the state of a physical system based on a set of noisy measurements is known as filtering. The state describes what a system is "doing" at any given time. The flight path of an aircraft, for example, may be represented by its position, velocity, and heading; position observations can be made in real time to estimate the velocity and heading of the aircraft. Bayesian state-space estimation methods, such as the Kalman filter and its variants, address the filtering problem with the aim of estimating the true state of a system. The adaptive approach predicts future states through a model of system equations and updates with respect to an observation. Links between state variables defined in model equations allow unobserved state values to also be updated.

This stochastic filtering approach suits the modeling of meteoroid trajectories using noisy and uncertain measurements. Typical meteoroid models mostly rely on measurements of the meteor/fireball brightness (Murray et al. 2000; Ceplecha & Revelle 2005; Kikwaya et al. 2011), though light curves tend to be variable and do not represent typical values predicted by single-body ablation models (Campbell-Brown & Koschny 2004). The meteoroid problem is complicated not only by unpredictable gross fragmentation in the atmosphere but also because the majority of initial state parameters are entirely unknown (m₀, σ, κ). Multiple approaches have been taken to handle these unknowns in fireball trajectory analysis. The manually intensive method of Revelle (2007) is based on the brute-force least-squares approach of Ceplecha & Revelle (2005). It does include the luminosity of the fireball (derived from manual interpretation of a light curve) as a proxy for mass loss and solves for fragmentation, as well as σ and κ. As it is still based on a least-squares optimization, model and observation errors are not rigorously examined; rather, overall errors are given as the standard deviation of residuals. The amount of manual input required also limits the number of fireballs that may be analyzed. The DFN observed over 300 fireball events in 2015 over its 2.5 million km² double station viewing area. This continental scale deployment of $\gt 50$ automated observatories has been possible owing to the low cost of each system. At this time, there is no expensive, high-voltage photomultiplier tube to measure fireball brightnesses. A trajectory analysis approach that is able to determine meteoroid parameters without a light curve, and which can be automated, will allow an unbiased evaluation of all events.

Very few models exist that enable the reduction of fireball data without a light curve. The method of Gritsevich (2009) solves for two dimensionless parameters rather than multiple unknown trajectory parameters. This still requires an initial accurate velocity and struggles with highly scattered data sets (Sansom et al. 2015). The various Kalman filtering methods used by Sansom et al. (2015, 2016) are fully automated techniques of determining the statistical likelihood of meteoroid state throughout bright flight and allow a robust analysis of observation and model errors. As with previous dynamical approaches to fireball modeling, these require a predetermined initial parameter set, withholding a general solution. To remove this limitation and fully analyze the statistical likelihood of the final state of a meteoroid given a range of likely initial states, we can use a method that combines a Monte Carlo (MC) approach with the filtering problem—a particle filter (Gordon et al. 1993). Simply, a "cloud" of particles are initiated with state values determined by a probability function. The "cloud" will be denser where probabilities are higher. Particles are propagated forward in time according to the state equations and weighted according to an observation. A new generation of particles are resampled from the existing pool, based on their weighting, and particles that are of low probability are preferentially removed.

The Bunburra Rockhole fireball was observed over the Australian outback by the DFN in 2007 and produced the network's first recovered meteorite (Spurný et al. 2012). An extended Kalman filter (EKF; Sansom et al. 2015) and an unscented Kalman filter (Sansom et al. 2016) have been used to model the Bunburra Rockhole fireball given a set of starting parameters. Neither filter explicitly includes gross fragmentation; however, Sansom et al. (2016) applied two unscented Kalman filters in an interactive multiple model (IMM) to determine likely periods of fragmentation. Here we will examine the suitability of this sequential MC technique for modeling fireball meteoroid trajectories using the Bunburra Rockhole fireball data set.

The technique used in this paper for estimating meteoroid parameters is one of a broader class of techniques known as Bayesian state-space methods. These methods involve encapsulating the knowledge of a system based on its state, given by the vector ${\boldsymbol{x}}$. The state of an object could be its position and velocity, for example. The probability of the object being in state ${\boldsymbol{x}}$ at time instant t_k is represented as the conditional probability density function

$$\begin{eqnarray}&&p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{z}}}_{1:k}),\end{eqnarray} \tag{ 1 }$$

where ${{\boldsymbol{z}}}_{k}$ is the observation of the system made at time t_k and ${{\boldsymbol{z}}}_{1:k}$ is the history of all observations up until time t_k.

The calculation of Equation (1) is achieved recursively through the application of Bayes's rule

$$\begin{eqnarray}&&p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{z}}}_{1:k})=\displaystyle \frac{p({{\boldsymbol{z}}}_{k}| {{\boldsymbol{x}}}_{k})p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{z}}}_{1:k-1})}{p({{\boldsymbol{z}}}_{k}| {{\boldsymbol{z}}}_{1:k-1})}.\end{eqnarray} \tag{ 2 }$$

The terms in the numerator of Equation (2) are defined through the state-space equations, while the denominator can simply be considered as a normalizing constant.

There are three state-space equations. The state prior initializes the recursion and encapsulates all prior information about the state of the system,

$$\begin{eqnarray}&&p({{\boldsymbol{x}}}_{0}).\end{eqnarray} \tag{ 3 }$$

The measurement equation relates the observations (e.g., position) to the state of the system (e.g., position and velocity)

$$\begin{eqnarray}&&{{\boldsymbol{z}}}_{k}=h({{\boldsymbol{x}}}_{k},{{\boldsymbol{w}}}_{k}),\end{eqnarray} \tag{ 4 }$$

where ${{\boldsymbol{w}}}_{k}$ is a stochastic noise process with known distribution. Equation (4) defines the likelihood function, $p({{\boldsymbol{z}}}_{k}| {{\boldsymbol{x}}}_{k})$, which is the first term in the numerator of Equation (2). The process equation models how the state evolves in discrete time,

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{k+1}=f({{\boldsymbol{x}}}_{k},{{\boldsymbol{u}}}_{k}),\end{eqnarray} \tag{ 5 }$$

where ${{\boldsymbol{u}}}_{k}$ is another noise process with known distribution. Equation (5) defines the transition density $p({{\boldsymbol{x}}}_{k+1}| {{\boldsymbol{x}}}_{k})$, which is incorporated into the second term in the numerator of Equation (2) through the Chapman–Kolmogorov equation (Jazwinski 1970),

$$\begin{eqnarray}&&p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{z}}}_{1:k-1})=\int p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{x}}}_{k-1})p({{\boldsymbol{x}}}_{k-1}| {{\boldsymbol{z}}}_{1:k-1})d{{\boldsymbol{x}}}_{k-1}.\end{eqnarray} \tag{ 6 }$$

This section outlines the state space and the state-space equations chosen to model the motion and measurement of a meteoroid process for the purposes of this paper. The specific parameters used in the model to estimate the trajectory characteristics of the Bunburra Rockhole data set are given in Section 5.

The state that defines the meteoroid system includes the physical parameters of motion, as well as trajectory parameters σ and κ:

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{k}=\left[\begin{array}{c}{l}_{k}\\ {v}_{k}\\ {m}_{k}\\ {\sigma }_{k}\\ {\kappa }_{k}\end{array}\right]\mathrm{position}\\ &&\mathrm{velocity}\\ &&\mathrm{mass}\\ &&\mathrm{ablation}\ \mathrm{parameter}\\ &&\mathrm{shape} \mbox{-} \mathrm{density}\ \mathrm{parameter},\end{eqnarray} \tag{ 7 }$$

where the position is measured along a predefined path produced by triangulating observations from several imaging sensors.

The measurement in Equation (4) is given by

$$\begin{eqnarray}&&{{\boldsymbol{z}}}_{k}=H{{\boldsymbol{x}}}_{k}+{{\boldsymbol{w}}}_{k},\end{eqnarray} \tag{ 8 }$$

where the measurement matrix is

$$\begin{eqnarray}H=\left[\begin{array}{ccccc}1 & 0 & 0 & 0 & 0\end{array}\right]\end{eqnarray} \tag{ 9 }$$

and the measurement noise process, ${{\boldsymbol{w}}}_{k}$, is Gaussian with zero mean and variance R_k.

As a meteoroid passes through the atmosphere, its behavior can be modeled by the aerodynamic equations from the single-body theory of meteoroid entry (Hoppe 1937; Baldwin & Sheaffer 1971) given by Equation (11), which uses atmospheric densities, ${\rho }_{a}$, acquired using the NRLMSISE-00 atmospheric model (Picone et al. 2002), local acceleration due to gravity, g, and entry angle from horizontal, ${\gamma }_{e}$. It is natural to model the change of meteoroid state as a continuous-time differential equation

$$\begin{eqnarray}&&\dot{{\boldsymbol{x}}}={f}_{c}({\boldsymbol{x}})+{{\boldsymbol{u}}}_{c},\end{eqnarray} \tag{ 10 }$$

where ${f}_{c}({\boldsymbol{x}})$ is defined using

$$\begin{eqnarray}&&\displaystyle \frac{{dl}}{{dt}}=v,\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}\displaystyle \frac{{dv}}{{dt}} & = & -\displaystyle \frac{1}{2}\displaystyle \frac{{c}_{d}A{\rho }_{a}}{{\rho }_{m}^{2/3}}{v}^{2}{m}^{(\mu -1)}+g\sin {\gamma }_{e}\\ & = & -\kappa {\rho }_{a}{v}^{2}\,{m}^{(\mu -1)}+g\sin {\gamma }_{e},\end{eqnarray} \tag{ 11b }$$

$$\begin{eqnarray}\displaystyle \frac{{dm}}{{dt}} & = & -\displaystyle \frac{1}{2}\displaystyle \frac{{c}_{h}A{\rho }_{a}}{{H}^{* }}{v}^{3}{m}^{\mu }\\ & = & -\sigma \kappa {\rho }_{a}{v}^{3}\,{m}^{\mu },\end{eqnarray} \tag{ 11c }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\sigma }{{dt}}=0,\end{eqnarray} \tag{ 11d }$$

$$\begin{eqnarray}&&\displaystyle \frac{d\kappa }{{dt}}=0,\end{eqnarray} \tag{ 11e }$$

and the continuous-time process noise, ${{\boldsymbol{u}}}_{c}$, is Gaussian with zero mean and covariance ${{\boldsymbol{Q}}}_{c}$. Time integration of Equation (10) is needed to arrive at the form required by the filtering state-space Equation (5). In this case

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{k+1}={\int }_{{t}_{k}}^{{t}_{k+1}}{f}_{c}({\boldsymbol{x}}){dt}+{{\boldsymbol{u}}}_{k}.\end{eqnarray} \tag{ 12 }$$

Due to the nonlinearities of Equation (11), the discrete-time process noise, ${{\boldsymbol{u}}}_{k}$, is not Gaussian, but can be closely approximated by Gaussian noise with zero mean and covariance

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{k}={\int }_{{t}_{k}}^{{t}_{k+1}}{e}^{{Ft}}\,{{\boldsymbol{Q}}}_{c}\,{e}^{{F}^{T}t}{dt}\end{eqnarray} \tag{ 13 }$$

(Grewal & Andrews 1993), where the matrix F is the linearized form of the process equation

$$\begin{eqnarray}&&F=\displaystyle \frac{\partial {f}_{c}({\boldsymbol{x}})}{\partial {\boldsymbol{x}}}.\end{eqnarray} \tag{ 14 }$$

Due to the form of the nonlinear functions given by Equation (11), the integrations required by Equations (12) and (13) cannot be found analytically. Numerical methods are used to calculate the integrals.

There are a range of methods for finding the distribution of ${{\boldsymbol{x}}}_{k}$ by solving Equation (2). The applicability of the method depends on the form of the state-space equations. If the measurement function and process function are linear and all the noise and prior distributions are Gaussian, then the solution to Equation (2) can be found analytically. This solution is known as the Kalman filter (Grewal & Andrews 1993). In the case where the equations are nonlinear or the distributions are non-Gaussian, such as the single-body equations for modeling meteoroid trajectory given by Equation (11), there are no exact solutions and approximations are required.

The EKF (Sansom et al. 2015) approximates the noise distributions as Gaussian and finds a linear approximation to the process equations. The unscented Kalman filter (Sansom et al. 2016) approximates the posterior distribution as a Gaussian, but avoids approximating the measurement or process equations through a method of statistical linearization (Särkkä 2007).

A particle filter does not require any assumptions about the form of the state equations or have any limitations on the noise distributions. This flexibility is achieved by representing the posterior density Equation (2) as a set of N_s weighted particles, which are simply points in the state space (Gordon et al. 1993; Arulampalam et al. 2002). The ith random particle at time t_k is represented by its state, ${{\boldsymbol{x}}}_{k}^{i}$, and weight, w_kⁱ,

$$\begin{eqnarray}&&\{{{\boldsymbol{x}}}_{k}^{i},{w}_{k}^{i}\}\quad i=1\,,...,{N}_{s}.\end{eqnarray} \tag{ 15 }$$

Weights are normalized so that

$$\begin{eqnarray}&&\displaystyle \sum _{i}^{{N}_{s}}{w}_{k}^{i}=1.\end{eqnarray} \tag{ 16 }$$

The probability distribution of the state is approximated by this set of weighted particles

$$\begin{eqnarray}&&p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{z}}}_{1:k})\approx \displaystyle \sum _{i}^{{N}_{s}}\delta ({{\boldsymbol{x}}}_{k}-{{\boldsymbol{x}}}_{k}^{i}){w}_{k}^{i},\end{eqnarray} \tag{ 17 }$$

where $\delta ({\boldsymbol{y}})$ is the Dirac delta function, defined such that

$$\begin{eqnarray}\delta ({\boldsymbol{y}})=\left\{\begin{array}{ll}1 & {\boldsymbol{y}}=0\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 18 }$$

Statistics can be computed on this set of particles, for example, the mean of the distribution at any time t_k is approximated by

$$\begin{eqnarray}&&{\hat{{\boldsymbol{x}}}}_{k}=\displaystyle \sum _{i}^{{N}_{s}}{w}_{k}^{i}{{\boldsymbol{x}}}_{k}^{i},\end{eqnarray} \tag{ 19 }$$

with the state covariance calculated as

$$\begin{eqnarray}&&\mathrm{Cov}({{\boldsymbol{x}}}_{k})=\displaystyle \sum _{i}^{{N}_{s}}{w}_{k}^{i}({{\boldsymbol{x}}}_{k}^{i}-{\hat{{\boldsymbol{x}}}}_{k}){({{\boldsymbol{x}}}_{k}^{i}-{\hat{{\boldsymbol{x}}}}_{k})}^{T}.\end{eqnarray} \tag{ 20 }$$

There are strong similarities between the implementation of a particle filter and the simpler Kalman filter. Both follow the following three steps:

• 1.  
  Initialization: start the filter with a known prior distribution, $p({{\boldsymbol{x}}}_{0}).$
• 2.  
  Prediction: propagate the distribution from time $k-1$ to time k using the process Equation (5).
• 3.  
  Update: use the measurement Equation (4) to update the predicted distribution with the measurement information, producing the posterior distribution at time k, $p({{\boldsymbol{x}}}_{k}| {{\boldsymbol{z}}}_{1:k}).$

The Kalman filter achieves these steps by exact analytic equations that manipulate the mean and covariance of the distribution at each step. On the other hand, the particle filter proceeds through calculation on each of the particles individually.

To initialize the particle filter, a set of particles are randomly sampled from the prior distribution, $p({{\boldsymbol{x}}}_{0})$, and weighted equally as ${w}_{0}^{i}=\tfrac{1}{{N}_{s}}$.

In the prediction step each particle is propagated forward in time via the process Equation (12). To incorporate the uncertainty of the system, a sample from the process noise, u_k, is randomly generated for each particle. Using the process equation to propagate the particles results in the simplest form of the filter. The particle filter literature generalizes this through importance sampling, where an arbitrary proposal distribution can be used, instead of the process equation (Arulampalam et al. 2002). Sophisticated proposal distributions can make a particle filter implementation more efficient (require fewer particles), but they have not been investigated for this application.

The update step adjusts the weight of each particle. The weight is obtained by evaluating the likelihood function for each particle

$$\begin{eqnarray}&&{\tilde{w}}_{k}^{i}=p({{\boldsymbol{z}}}_{k}| {{\boldsymbol{x}}}_{k}^{i}){w}_{k-1}^{i}.\end{eqnarray} \tag{ 21 }$$

The weights are then normalized to satisfy Equation (16),

$$\begin{eqnarray}&&{w}_{k}^{i}=\displaystyle \frac{{\tilde{w}}_{k}^{i}}{{\displaystyle \sum }_{i}^{{N}_{s}}{\tilde{w}}_{k}^{i}}.\end{eqnarray} \tag{ 22 }$$

Over time the particle weights can transfer to a few select particles, thereby updating insignificant particles at the expense of computing power (Arulampalam et al. 2002). This is known as the degeneracy problem, and Equation (23) gives an approximate measure of particle effectiveness that can be used to assess the severity of the issue (Arulampalam et al. 2002):

$$\begin{eqnarray}&&{\hat{N}}_{\mathrm{eff}}={\left(\displaystyle \sum _{i}^{{N}_{s}}{({w}_{k}^{i})}^{2}\right)}^{-1}.\end{eqnarray} \tag{ 23 }$$

The degeneracy problem can be addressed by resampling the data after weights have been calculated. A new population of particles is generated from the current sample pool based on given weightings, the objective being to preferentially remove samples of lower weights. The probability of resampling any given particle i is w_kⁱ. The optional resampling step is taken if the number of effective particles drops below some threshold. After resampling, all of the particle weights are set to $1/{N}_{s}$.

Dedicated fireball networks, such as the DFN, capture fireball events from multiple locations, providing triangulated position observations with time. This also enables a rough calculation of velocities throughout the trajectory.

5.1. Initialization

When initializing the state prior for the set of N_s particles at the start of the luminous trajectory (t₀), the initial position and, to an extent, the initial velocity⁴ can be reasonably well constrained. The other state parameters, $m,\sigma ,\kappa$, however, are not directly observable. To explore the data space and determine likely values for m₀, as well as constants σ and κ, each particle is initiated with a random value within a given range. The state prior for each particle is initialized according to Table 1, with m₀^min in all cases set to 0.5 kg.

Table 1.  Description of the Method Used by the Particle Filter to Initialize State Parameters for Each Particle

Parameter to be Initiated	Method Used
l0	Random choice based on Gaussian ${ \mathcal N }(0,10\,{\rm{m}})$
	(from triangulation errors)
v0	Random choice based on Gaussian ${ \mathcal N }({v}_{0},500\,{\rm{m}}\,{{\rm{s}}}^{-1})$
	(from triangulation errors)
m0	Random choice from 0 to m0max (kg)
σ	Random choice between 0.001 and 0.05 ${{\rm{s}}}^{2}\,{\mathrm{km}}^{-2}$
	(from Ceplecha et al. 1998 for asteroidal material)
κ	cd—random choice based on Gaussian ${ \mathcal N }(1.3,0.3)$
	(based on aerodynamic drag values from Zhdan et al. 2007)
	A—random choice based on Gaussian ${ \mathcal N }(1.4,0.33)$
	(close to spherical values)
	${\rho }_{m}$—the PDF representing meteorite bulk densities is multimodal; to fully represent this distribution, initialization is performed in two stages
	First, a random choice of meteorite type is made based on recovered percentages (80% chondrites, 11% achondrites, 2% stony-iron, 5% iron, 2% cometary; Grady 2000)
	Second, a random choice of bulk density is made based on the Gaussian PDF representing chosen meteorite type;
	chondrites—${ \mathcal N }(2700,420)$ (after Britt & Consolmagno 2003);
	achondrites—${ \mathcal N }(3100,133)$ (after Britt & Consolmagno 2003);
	stony-iron—${ \mathcal N }(4500,133)$ (after Britt & Consolmagno 2003);
	iron—${ \mathcal N }(7500,167)$ (after Consolmagno & D. T. Britt 1998);
	cometary—${ \mathcal N }(850,117)$ (after Weissman & Lowry 2008).

Note. A random selection is made for each value using either a Gaussian probability density function (PDF) (mean and standard deviation given), a uniform PDF within a given value range, or a multimodal distribution in the case of bulk density.

Download table as:  ASCIITypeset image

5.2. Prediction

At every observation time, t_k, the state of each particle is evaluated using the system model given by Equation (10). ${{\boldsymbol{Q}}}_{c}$ values used here to represent the continuous process noise in the given model for meteoroid trajectories are given by the diagonal matrix represented by Equation (24). The diagonal elements of ${{\boldsymbol{Q}}}_{c}$ in Equation (24) are the variance values for dl/dt, dv/dt, dm/dt, $d\sigma /{dt}$, and $d\kappa /{dt}$, respectively. The uncertainties in position and velocity are introduced through noise in the acceleration model given by Equation (11(b)), and the variance for dl/dt for this process model is therefore set to $0\,{\rm{m}}\,{{\rm{s}}}^{-1}$. The other model equations, however, are not able to represent the system in its entirety; complications, such as fragmentation, affect all other state process models. At this stage, we assume that the shape density and ablation parameters will not change dramatically over the meteoroid flight and are attributed to small process noise values. There is a high uncertainty in the mass loss for the single-body ablation model given by Equation (11(c)), and so a large range of masses are allowed to be explored by the particles. The process noise in mass is a multiple of the mass in order to keep it within a consistent order of magnitude. The discrete process noise, ${{\boldsymbol{Q}}}_{k}$, is calculated at every time step following Equation (13):

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{c}=\text{diag}{\displaystyle \left[\phantom{\tfrac{1}{2}}(0\text{m}{\text{s}}^{-1}),(75\text{m}{\text{s}}^{-2}),(0.2\times {m}_{k}\text{kg}{\text{s}}^{-1}),({10}^{-4}\text{s}{\text{km}}^{-2}),({10}^{-5}{\text{m}}^{2}{\text{kg}}^{-2/3}{\text{s}}^{-1})\right]}^{2}\end{eqnarray} \tag{ 24 }$$

To improve the computation time of this method, the nonlinear integration (12) of all N_s particles, as well as their associated ${{\boldsymbol{Q}}}_{k}$, is performed simultaneously using parallel multiprocessing.

5.3. Update

The triangulated position of the meteoroid along the trajectory at time k is the observation measurement ${{\boldsymbol{z}}}_{k}$. The weight $({\tilde{w}}_{k}^{i})$ for each particle ${{\boldsymbol{x}}}_{k}^{i}$ is calculated using a 1D Gaussian probability distribution function

$$\begin{eqnarray}&&p({{\boldsymbol{z}}}_{k}| {{\boldsymbol{x}}}_{k}^{i})=\displaystyle \frac{1}{\sqrt{2{R}_{k}\pi }}{e}^{-\tfrac{{({{\boldsymbol{z}}}_{k}-{{\boldsymbol{Hx}}}_{k}^{i})}^{2}}{2{R}_{k}}}\end{eqnarray} \tag{ 25 }$$

in Equation (21), with the observation noise having a variance ${R}_{k}={(100{\rm{m}})}^{2}$. This is based on errors in timing and triangulated position, reflecting the accuracy of the data set being used.

In order to avoid degeneracy in the particle set, we have use the stratified resampling method described by Arulampalam et al. (2002) after each update step.

The data acquired by Spurný et al. (2012) for the Bunburra Rockhole fireball are used to test the suitability of the particle filter in estimating the state of a meteoroid during atmospheric entry. The Bunburra Rockhole data set consists of 113 published observations of position with time along the trajectory. Note that no observation data were published between t = 0.0 s and t = 0.1899 s or from t = 5.3165 s to t = 5.4589 s. Our modeling will use times relative to t₀ = 0.1899 s along the trajectory. A particle filter is run using a set of 10,000 particles (${N}_{s}={\rm{10,000}}$). Particles are initiated according to Table 1, with m₀^max set to 2000 kg.

Figure 1 shows all the resulting particle masses with weights $\gt 0$ from t₀ to t_end. The range of σ and κ values used to initiate each particle results in a variety of predicted trajectory "paths."

Zoom In Zoom Out Reset image size

To aid in understanding the different trajectories predicted by the particle filter, five particles at t₀ have been selected for discussion (${{\boldsymbol{x}}}_{0}^{j}$ given in Table 2). Figure 2 highlights these particles, ${{\boldsymbol{x}}}_{0}^{a-e}$, along with all particles that are generated from them at later time steps (by either propagation from ${t}_{k-1}$ or resampling at t_k).

Zoom In Zoom Out Reset image size

Table 2.  State of Five Particles at t₀

${{\boldsymbol{x}}}_{0}^{j}$	l0	v0	m0	${\sigma }_{0}$	${\kappa }_{0}$	Reference Color
	(m)	$(\mathrm{km}\,{{\rm{s}}}^{-1})$	(kg)	$({{\rm{s}}}^{2}\,{\mathrm{km}}^{-2})$	(m2 kg−2/3)	in Figure 2
${{\boldsymbol{x}}}_{0}^{a}$	−1.57	12.80	10.1	0.022	0.0083	blue
${{\boldsymbol{x}}}_{0}^{b}$	−18.60	12.88	14.3	0.020	0.0058	green
${{\boldsymbol{x}}}_{0}^{c}$	5.00	12.48	176.2	0.021	0.0039	red
${{\boldsymbol{x}}}_{0}^{d}$	−17.19	12.96	212.1	0.037	0.0083	dark orange
${{\boldsymbol{x}}}_{0}^{e}$	12.41	13.10	234.0	0.041	0.0133	light orange

Note. All future particles resampled from these are highlighted in Figure 2 according to the color given here.

Download table as:  ASCIITypeset image

The variation in σ (Figure 2(b)) and κ (Figure 2(c)) values with time is due to the addition of process noise, ${{\boldsymbol{u}}}_{k}$, in Equation (10). As this noise is random Gaussian, it allows small variations between identical resampled particles that would have originally shared equal values. Areas of greater particle density are characteristic of higher probability states.

Orange particles in Figure 2 originate from ${{\boldsymbol{x}}}_{0}^{e}$. The steep change in mass with time (Figure 2(a) is due to the high σ (Figure 2(b)) and κ (Figure 2(c) values with which they were initiated. Particles that no longer fit the observed data are preferentially removed by the resampling process, and their "path" discontinues in Figure 2. Although particles originating from ${{\boldsymbol{x}}}_{0}^{c-e}$ were initiated with diverse σ (Figure 2(b)) and κ (Figure 2(c)) values, they, along with all other particles with ${m}_{0}^{i}\gt 27\,\mathrm{kg}$, have insignificant weight past 5.0 s. A visual comparison of predicted particle velocities with velocities calculated from position measurements is shown in Figure 2(d). The "survival" of ${{\boldsymbol{x}}}_{0}^{a,b}$ to t_end is due to their higher w_kⁱ values, indicating superior fits to the observations (and visually noticeable in Figure 2(d)).

The final trajectory parameters of the Bunburra Rockhole meteoroid have been previously determined by Spurný et al. (2012) using the dynamic gross fragmentation model (GFM) of Ceplecha et al. (1993) and the meteoroid fragmentation model (MFM) of Ceplecha & Revelle (2005), which integrates fireball brightness with the dynamics (Table 3). Both the GFM and MFM require initial assumptions, including the entry mass and a manually predefined fragmentation pattern based on the light curve (Ceplecha & Revelle 2005). Errors given by these models relate to the standard deviation of the residuals between modeled and measured observations; observational uncertainties, assumptions made in the model, and model parameters are not propagated. The Kalman filter methods applied by Sansom et al. (2015, 2016) to meteoroid trajectory modeling perform a comprehensive analysis of the errors of both model and observations but share the limitations of previous models in requiring a single set of initial entry parameters to be predetermined.

Table 3.  Mean Final State Values Estimated by the Particle Filter (19), Alongside Published Values

Method	lend	vend	mend	${\sigma }_{\mathrm{end}}$	${\kappa }_{\mathrm{end}}$
	(km)	($\mathrm{km}\,{{\rm{s}}}^{-1}$)	(kg)	(${{\rm{s}}}^{2}\,{\mathrm{km}}^{-2}$)	(m2 kg−2/3)
GFMa			1.5 ± 0.2	0.0331 ± 0.0007
				(apparent)
MFMa		5.77 ± 0.04	1.1	$0.002\pm 0.001/0.004$	0.0035
				(intrinsic)
Dynamicb	60.07	6.109	2.36	0.0371	0.0062
optimization				(apparent)
EKFb	60.03 ± 0.062	6.05 ± 0.24	2.30 ± 1.63
UKFc	60.04 ± 0.058	6.10 ± 0.20	2.88 ± 1.04
IMMc	60.01 ± 0.007	5.90 ± 0.06	1.32 ± 0.49
Particle filter	59.89 ± 0.038	6.03 ± 0.22	2.16 ± 1.33	0.0219 ± 0.0007	0.0042
				(apparent)	± 0.000

Notes. Errors given by all previous methods reflect only model errors within the given initial input assumptions given. The GFM and MFM methods do not consider observation uncertainties (Ceplecha & Revelle 2005). The particle filter errors are calculated as $\sqrt{\mathrm{Var}({\hat{{\boldsymbol{x}}}}_{k})}$ given by Equation (20) and alone give a fully inclusive analysis of the trajectory model and observation uncertainties to provide a more realistic understanding of real-world variability.

^aSpurný et al. (2012); GFM = gross fragmentation model; MFM = meteoroid fragmentation model. ^bSansom et al. (2015); κ value determined using ${c}_{d}\,=\,1.3;$ EKF = extended Kalman filter. ^cSansom et al. (2016); UKF = unscented Kalman filter; IMM = interactive multiple model.

Download table as:  ASCIITypeset image

The statistical approach of the particle filter is not limited to any one set of input parameters. It encapsulates all prior knowledge of the parameter space by exploring the full range of plausible parameter values to produce an unbiased analysis. Given that model and observation uncertainties are incorporated and propagated, this method provides a statistically robust final state estimate that is no longer dependent on any single set of assumed input parameters, providing a more realistic understanding of real-world variability. The independence of the particle filter and lack of manual input enable full automation of this method.

The spread of final particle states at t_end can be summarized by the weighted mean (19) in Table 3. Errors are calculated as the square root of the covariance diagonal elements given by Equation (20). The ablation parameter is an interesting result. Although the particle filter does not explicitly model fragmentation, Q_c allows for a certain amount of variation in state parameters due to unmodeled processes and inherently includes fragmentation to some extent, without the need for a predefined fragmentation pattern (required by MFM; Ceplecha & Revelle 2005). As discussed by Ceplecha & Revelle (2005), the intrinsic value of the ablation parameter remains constant throughout the trajectory regardless of fragmentation. When fragmentation is not modeled explicitly, variations in the ablation parameter appear to occur and must therefore be expressed as the apparent ablation parameter. The GFM produces an apparent σ, whereas the MFM, as it incorporates the light curve, is able to define the intrinsic σ. The value determined using the particle filter is slightly lower than the apparent σ of the GFM, and it is therefore plausible that we can use this difference to quantify the extent to which fragmentation is included in the final state estimate.

Using a particle filter, the state estimates at each time step are iteratively updated based on the past data; future observations are not included. The final states alone result from processing all observations. As a predicted particle becomes inconsistent with the observations, it becomes an unlikely scenario for future times, but it does not mean that this original path can be discounted. It is noticeable at ${t}_{k}=3.32\,{\rm{s}}$ that there is a drastic reduction in the number of particle "paths" that fit the observational data. The parameter space after this time is much more constrained. All particles at t_end originate from particles with ${{\boldsymbol{x}}}_{0}\lt 27\,\mathrm{kg};$ these particles are consistent with both parts of the trajectory displaying no dramatic change in mass. It is possible that particles of initially higher mass are discontinued in favor of lower-mass scenarios as a result of gross fragmentation reflected in the observation data. Without including all the data at every time step, the most likely state 'path' for the entire trajectory cannot be constrained; we cannot distinguish the full particle history.

In order to distinguish likely initial masses, we need to be able to explore drastic changes in mass. The IMM smoother as described by Sansom et al. (2016) has this capability and uses all observational data at each time step. It, however, requires a single predefined set of initial parameters. This is a well-suited complementary method to our current implementation of a particle filter. The particle filter framework, however, is flexible enough to incorporate dynamic models that explicitly capture gross fragmentation events. Future work will explore more sophisticated dynamic models, as well as particle filter smoothing, to reconstruct the full meteoroid trajectory.

Including brightness as a state in trajectory modeling would also provide an additional observation with which to weight particles. As brightness is linked to mass, its addition would not only improve state estimates but also inherently include information on fragmentation.

The use of a particle filter to approximate fireball trajectories provides a statistical analysis of the meteoroid state, including unobservable trajectory parameters. This is the first approach of its kind in this field. Other nonlinear filtering algorithms such as the EKF (Sansom et al. 2015) and the unscented Kalman filter (Sansom et al. 2016), as well as other least-squares approaches (Ceplecha et al. 1993; Ceplecha & Revelle 2005), require a predetermined set of initial parameters to statistically analyze the trajectory of a meteoroid. The iterative MC simulations of a particle filter not only are capable of automating the analysis of fireball trajectories but also are able to do so without the need for limiting input parameters to single assumed values; rather, it encapsulates all prior knowledge of the parameter space, to produce an unbiased analysis. The adaptive filter approach uses the observations of the meteoroid's position as it travels through Earth's atmosphere to update state estimates. Predicted positions similar to those observed are given a higher weighting and are preferentially resampled at the next time step. This gives a final state estimate (Table 3) with robust error propagation of uncertainties in the initial parameters, observations, and the dynamic model (e.g., unpredictable gross fragmentation events). Even though trajectory parameters σ and κ are not currently set to vary systematically with time (noise is added to create diversity between resampled particles to avoid degeneracy only), a stochastic approach to their determination has not previously been conducted. Incorporating brightness as an additional state will provide supplementary data and improve estimates. This method currently allows an automated dynamic analysis of fireball trajectories.

This work was funded by the Australian Research Council as part of the Australian Laureate Fellowship scheme.