Inference of meteoroid characteristics using a genetic algorithm

Icarus 329 (2019) 270–281

Contents lists available at ScienceDirect

Icarus journal homepage: www.elsevier.com/locate/icarus

Inference of meteoroid characteristics using a genetic algorithm a,d,⁎

Ana María Tárano

c

d

, Lorien F. Wheeler , Sigrid Close , Donovan L. Mathias

b

T

a

Science and Technology Corporation, NASA Ames Research Center, MS 258, Moffett Field, CA 94035, United States NASA Ames Research Center, MS 258-5, Moffett Field, CA 94035, United States RedLine Performance Solutions, NASA Ames Research Center, MS 258-6, Moffett Field, CA 94035, United States d Stanford University, Stanford, CA 94305, USA b c

ARTICLE INFO

ABSTRACT

Keywords: Asteroid Genetic algorithm Meteoroid Impact risk Asteroid characterization

A methodology is introduced to optimize and extend the inference of pre-entry size, density, strength, and mass of asteroids based on observed light curves. In this development study, a genetic algorithm (GA) approach is coupled with the fragment-cloud model (FCM) to efficiently evaluate entry and breakup for numerous potential asteroid property combinations and determine which case best matches the observed data. FCM produces energy deposition curves based on assumed pre-entry conditions, and the GA finds values for these inputs that minimize an objective function characterizing the difference between the FCM curve and a target curve. We present an overview of the GA approach, and then demonstrate its capability to infer pre-entry properties for three wellcharacterized events: Chelyabinsk, Lost City, and Benešov. In all cases, our initial mass and size estimates were within the range of published values.

1. Introduction Quantifying the asteroid impact risk to the Earth requires specific knowledge about the object's pre-entry properties. In the event of a known, specific impactor, it is possible to determine detailed astronomical measurements of the trajectory, brightness, and diameter using thermal and optical observations. However, it is unlikely that such measurements will be available early in the threat detection process, especially for small impactors. To calculate the ensemble risk from the overall asteroid population or an impactor with uncertain properties, impact risk assessment studies such as Mathias et al. (2017) and Rumpf et al. (2017) use statistical distributions to describe the range and likelihood of potential physical properties. But even in the development of those distributions, obvious difficulties exist due to the lack of detailed direct measurements (Borovička et al., 2017; Brown et al., 2016; Nemtchinov et al., 1997; Pravec et al., 2012). As a result, the distributions must be formed by combining various sources of available information. Such sources can include astronomical measurements, rendezvous spacecraft mission data, meteorite measurements, or measurements of bolides as they pass through the atmosphere. Ceplecha et al. (1998) present an excellent overview of the historical use of instruments to study meteors. Bolide measurements, typically in the form of light curves, are becoming more readily available through the presence of space assets (Tagliaferri et al., 1994) and growing ground-

⁎

based camera networks (e.g., Cooke and Moser, 2012). These measured light curves provide a record of the energy that the bolide deposits in the atmosphere as it fragments, ablates, and decelerates during entry. This data, along with any available trajectory or meteorite data, can be used to study breakup characteristics, refine entry modeling approaches or parameters, and make inferences about key pre-entry properties such as size, mass, density, and strength. A variety of meteor modeling studies have endeavored to infer or bound unknown properties by using theoretical, semi-analytic, or numerical models of entry and breakup physics to reproduce observed light curves, the associated atmospheric energy deposition, and/or dynamic trajectory data (e.g., Borovička et al., 1998; Ceplecha and ReVelle, 2005). The suitability of the potential initial properties or modeling assumptions are generally evaluated by matching modeled fragmentation events to peaks in the observed light curve through trialand-error fits (e.g., Borovička et al., 1998, 2013, 2015, 2017; Gi, 2017). These studies have employed a number of different semi-analytic fragmentation modeling approaches, including single-body “pancake” or “liquid drop” models (e.g., Chyba et al., 1993; Hills and Goda, 1993; Melosh, 1981; Zahnle, 1992), progressive fragmentation approaches (e.g., Baldwin and Sheaffer, 1971; Ceplecha and ReVelle, 2005; Levin, 1956; ReVelle, 2007), and hybrid models that combine multiple approaches (e.g., Wheeler et al. (2017, 2018); Popova et al. (2013), Supplement; Borovička et al., 2013). Register et al. (2017) provide a

Corresponding author at: Science and Technology Corporation, NASA Ames Research Center, MS 258, Moffett Field, CA 94035, United States. E-mail address: [email protected] (A.M. Tárano).

https://doi.org/10.1016/j.icarus.2019.04.002 Received 24 August 2018; Received in revised form 20 February 2019; Accepted 3 April 2019 Available online 10 April 2019 0019-1035/ © 2019 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/).

Icarus 329 (2019) 270–281

A.M. Tárano, et al.

comparison of several of the predominant approaches. The results of such models have been used to estimate general properties such as total mass, energy, or size (e.g., Brown et al., 2002, 2013); to evaluate aspects of the entry physics such as ablation rates and luminous efficiency (e.g., Artemieva and Shuvalov, 2001; Ceplecha and ReVelle, 2005; ReVelle, 2007); to assess breakup characteristics and aerodynamic strengths (e.g., Borovička et al., 2013, 2017; Brown et al., 2013); to reconstruct specific fragmentation sequences and predict masses and locations of meteorite falls (e.g., Artemieva and Shuvalov, 2001; ReVelle, 2007); and to infer structural characteristics (e.g., Borovička et al., 2017; Wheeler et al., 2018). However, the model matching process typically used in such studies is arduous, requiring extensive manual trial-and-error variations to be assessed. As a result, it can be difficult to fully explore the broad range of potential properties and modeling parameters, and the majority of observed meteors are not modeled unless meteorites are found and trajectory is observed. Furthermore, the matches obtained can be extremely dependent on case-specific modeling assumptions (e.g., Ceplecha and ReVelle, 2005), solutions can be non-unique (e.g., Popova et al. (2013), Supplement), and the relative quality of various model fits can be subjective or difficult to compare. Many authors—such as Borovička et al. (2013; Borovička et al. (2017); Brown et al. (2002, 2013); Popova et al. (2013); Wheeler et al. (2017); Wheeler et al. (2018)—allude to the ambiguity in their modeling matches because different combinations of properties or model parameters can lead to similar light curve or energy deposition profiles. Uncertain modeling parameters and varying degrees of model complexity pose additional barriers and can amplify these challenges. Many meteor modeling parameters (such as ablation coefficients, spread rates, or luminous efficiency) remain notably uncertain and poorly characterized due to lack of definitive data, and appropriate values may vary widely depending on both the modeling approach and the specific observed meteor. Similarly, more complex models can better match a realistic range of light curve features, potentially supporting more refined inferences, but also introduce more free parameters into the problem. Together, these challenges have limited the potential inference that could be gleaned from the growing availability of light curve data. In order to optimize and expand future inference capabilities, there is a need for a systematic way of evaluating and comparing a large range of entry cases. In the present work, we endeavor to surmount some of the barriers of characterizing asteroid properties by automating the process of matching entry modeling results to observed light curve data in a comprehensive and systematic manner. An automated tool for meteoroid entry modeling addresses the challenges mentioned above by enabling a thorough exploration and evaluation of potential solutions and uncertain modeling parameters to more effectively infer initial meteoroid properties. We have developed such a capability by combining a genetic algorithm (GA) approach with a semi-analytic asteroid entry and breakup model, known as the fragment-cloud model (FCM) (Wheeler et al., 2017). FCM provides the ability to represent a variety of breakup behaviors and light curve features without manually prescribing specifics about the break sequence (Wheeler et al., 2017, 2018), making it a good candidate for automated matching. The GA approach uses an evolution-based search algorithm, which finds solutions using principles of genetics and natural selection (Holland, 1984). For our meteor modeling application, the initial asteroid properties and modeling parameter values comprise the genes of the entry cases being evaluated and evolved. FCM is used to model each case in the population, and an objective function evaluates the fitness of the cases by comparing the FCM energy deposition results to those derived from an observed light curve. Subsequently, the GA evolves the selection of input parameters to find an optimal match. Together, the GA and FCM pair accomplish an automated search by comparing a comprehensive range of cases and strategically minimizing the error between the modeled and observed asteroid entry signatures.

In this paper, we present the current development of our GA design and the objective function best suited for matching observed meteor light curve data. We then demonstrate the approach's capability to automatically match a diverse range of curves and infer key parameters for three well-documented observed meteors. Results are presented for the 1970 Lost City fireball, the 1991 Benešov bolide, and the popularly recorded 2013 Chelyabinsk superbolide. With these examples, we aim to evaluate the most readily inferred parameters such as pre-entry mass, size, density, and strength. Lastly, we assess the advantages and challenges of the GA approach and define future work. 2. Methodology By combining the fragment-cloud model (FCM) and a genetic algorithm (GA), we enable effective automation of meteoroid energy deposition (light curve) matching by searching the parameter space in a systematic manner. Once an optimal match is obtained, the FCM input parameters that produced the optimal solution are used to infer the preentry state of the meteor. In this section, we first review FCM fundamentals and list the underlying assumptions used in the current approach. We then motivate the use of the GA and describe its structure, including the population characteristics, the genetic evolution process, and the objective function used to evaluate the fits. 2.1. Fragment-cloud model (FCM) The fragment-cloud model (FCM), described in detail in (Wheeler et al., 2017, 2018), is a semi-analytic model that simulates an asteroid's trajectory and fragmentation in order to compute the energy deposited in the atmosphere as a function of altitude. The model integrates the single-body equations of motion and ablation, and represents the breakup process using a combination of individual fragments and debris clouds. Fig. 1 shows a pictorial representation of the model with accompanying equations. The equations of motion depend on velocity v, diameter (included within cross-section A), mass m, flight angle θ, and atmospheric density ρair. The model assumes that the original body and fragments are spherical in shape. Mass loss depends on the ablation coefficient σ, atmospheric density, velocity, and diameter. A fragmentation event occurs whenever the stagnation pressure (ρairv2) exceeds the body's aerodynamic strength parameter S. This strength does not specifically represent tensile, compressive, or yield strength but rather is a proxy based on the aerodynamic load leading to the breakup. Upon each fragmentation, the body is split into a given number of child fragments and a debris cloud. The strengths of the child fragments S2 increase according to a Weibull-like size-strength scaling relation, controlled by the fragment strength scaling parameter α (Tsvetkov and Skripnik, 1991). The modeling parameters used to control the breakup behavior are the initial meteoroid strength, the strength-scaling coefficient, the number of fragments produced in each break, and the fraction of mass going into a debris cloud with each break, which we call cloud fraction. These fragmentation parameters, along with the ablation coefficient, are set for each case and remain constant throughout the entry. We also assume that masses of fragments are divided evenly at each break and both the cloud and the fragments have the same ablation coefficients. (Wheeler et al., 2017) provides a complete description of FCM and describes the model parameters in detail. 2.2. Genetic algorithm background Genetic algorithms are optimization tools that find a suitable solution by performing a directed search based on combinatorics (Blum and Roli, 2003). Like many optimization routines, the goal is to efficiently determine the best values for a set of decision variables that minimize a given evaluation function. The decision variables are permitted to take on a different combination of values from within the constraints of an 271

Icarus 329 (2019) 270–281

A.M. Tárano, et al.

Fig. 1. This image illustrates the parameters used in the fragment-cloud model (FCM) to simulate the entry and breakup of meteoroids entering the atmosphere. At each break, the body splits into a given number of child fragments and a debris cloud of a given mass fraction. The equations listed on the right side show how FCM computes the energy deposited by the meteoroid, with the air density ρair computed from a curve fit of the 1976 standard atmosphere tables (Wheeler et al., 2017).

A GA's ability to perform both global and local searches (i.e., explore the entire search space and refine solutions, respectively) is particularly useful for our problem because meteor modeling provides a non-linear and high-dimensional solution space (Mitchell, 1996). Furthermore, GAs are inherently appropriate for parallel computing, which dramatically decreases the turnaround time of the optimization process (Mitchell, 1996). A main objective of our work is to lower the time and effort of meteoroid modeling, which requires efficient computation times. We use GAlib as the tool for automating the meteoroid modeling process. The details behind this C++ library are described in Wall (1996). GAlib is well-documented and has been successfully used for many diverse applications across multiple industries (e.g., Demiriz et al., 1999; Huuk et al., 2014; Seppelt and Voinov, 2002). We minimally customized the GAlib to allow for parallel computing using Message-Passing Interface standards (Message Passing Interface Forum, 2015). Otherwise, we used the built-in GAlib features, which provided a rich set of GA implementations and their associated tuning parameters.

allowable domain. The function that the optimization program attempts to minimize is called the objective function. The set of all possible assignments of the objective function, determined by the exhaustive set of all combinations of decision variables that satisfy the constraints, is termed the solution space, search space, or feasible region. The GA performs an iterative search within this space to locate the best set of decision variables. GAs are part of a class of optimization known as evolutionary computation, where the search for the best solution is powered by principles of Darwinian natural selection and genetics. This approach enables the use of randomness in an intelligent and adaptive manner to minimize computational costs while converging to good—if not the best—solutions (Stützle, 1999). The basic structure of a GA is to initialize a set of potential solutions, known as a population, evaluate their fitness using the objective function, and evolve the population until a stopping criterion is met. This computational model of evolutionary processes includes abstractions of selection, mating, and mutation, which are defined by a series of genetic operators. These evolutionary forces drive the inspection of the search space by providing a bias to the kinds of solutions the optimization finds. A set of decision variables within a population is known as an individual, genome, or chromosome within the GA context (Holland, 1984). Each individual is composed of a set of genes, with each gene representing a particular decision variable within the individual. The selected value of a gene is also known as the allele of that gene. The genetic representation, also known as genetic encoding, corresponds to the design and implementation of the decision variables within a genome. The range of solutions that the GA can find depends on the genetic representation of the individuals because this defines the bounds of the search space. The objective function is evaluated to determine the fitness level of each individual within a population (Mitchell, 1996). This fitness level determines how likely an individual is to be selected to mate and sometimes mutate to make part of the next iterative generation. The stopping criterion is set by the GA designer. In this paper, we terminate the evolution process after a specified number of generations. When the evolution terminates, the best solution is returned to the user. Furthermore, one can save all of the populations' individuals, including their fitness level, to gain a better understanding of the GA's performance and evaluate what parameters are most readily inferred and which remain largely unconstrained.

2.3. Development of GA for meteor modeling We have adapted the genetic algorithm (GA) approach to solve the inverse problem of automatically inferring meteoroid characteristics from energy deposition data derived from optical light curves. In our meteor modeling application, the decision variables are the asteroid properties and fragmentation parameters that are used as inputs to model the entry in FCM. The individuals within a population are potential meteor entry cases, with each input parameter representing a gene on their chromosome. The objective function evaluates the fitness of the cases by computing the level of error between the resulting FCM energy deposition curve and that of the observed meteor. Genetic operators control the fitness-based selection, mating, and mutation of the cases to produce subsequent generations of cases for evaluation. The genetic makeup of the individuals that compose the population define the bounds of the search space while the genetic operators drive the inspection of the domain. Therefore, the design of both the encoding of the genes and the hyperparameters guiding the evolution (through the genetic operators) are responsible for the GA's success at approximating the global optimum. The GA hyperparameters included generation type and size, mutation and crossover approaches and probabilities, population size, and selection mechanism. 272

Icarus 329 (2019) 270–281

A.M. Tárano, et al.

including smaller meteoroids is beyond the scope of this study. Furthermore, for altitudes above 100 km, the ram pressure on the leading surface of a meteoroid is too small to disrupt even the most fragile of particles (Campbell-Brown and Koschny, 2004) so fragmentation is unlikely initialized above 100 km for larger meteoroids of interest. According to Silber et al. (2018), other flow regimes are not as well understood and smaller bodies cannot be treated like larger events.

We initially developed the GA methodologies and hyperparameters for our meteor modeling application using controlled experiments that tested their ability to match a set of synthetic FCM-generated curves with known parameters. This approach enabled us to prototype the GA independent from FCM modeling assumptions and challenges. The main objective of this prototyping was to determine the GA hyperparameters best-suited to evaluating meteoroid energy deposition curves. We manually tuned the GA hyperparameters until the randomly generated synthetic curves were robustly and quickly replicated, regardless of the initialized population. We used these curves to vet and design the GA hyperparameters by gradually increasing the number of FCM modeling parameters, in different combinations, that the GA would solve for. The synthetic curves demonstrated that energy tradeoffs could be mitigated by assuming entry trajectory information. The developed approach used to solve the synthetic curves is then applied to the observed meteor cases in the Results section. The choice of the GA's genetic operators, objective function, genetic representation, and population model is extremely dependent on the particular problem it attempts to solve (Gen and Cheng, 1997). Therefore, we used an FCM-generated set of synthetic curves with known inputs to establish the GA hyperparameters appropriate for our meteor modeling application. We randomly selected 5 sets of FCM initial conditions and generated energy deposition curves for each. We then manually varied the GA model hyperparameters until a robust set of matches was derived. The overall quality of the solutions was determined by considering the minimum values of the objective functions, qualitative assessments of the curve matches, and the ability to reproduce the known inputs with fastest computation times.

2.3.2. Fitness evaluation The objective function serves as the quantitative metric that represents the fitness of a given solution. In this case, the objective function is a measure of the difference between the FCM energy deposition curve and that of the bolide. Many different objective functions were tested to gain a better understanding of the limitations and strengths of different options, especially their behavior with regards to capturing bursts and flares. The tested objective functions are summarized below, and included comparisons of the energy deposition rate at each altitude, the cumulative energy deposited, and statistical associations. In the first type of function considered, we define the error εi at a given altitude i as

Bulk density (g/cm ) Initial breakup strength (kPa) Diameter (m) Ablation coefficient (kg/J) Cloud fraction Number of fragments Strength scaling

Minimum value

Maximum value

1.5 1 0.1 1e-9 0.1 2 0.1

5 10,000 100 9e-8 0.9 16 0.9

Ei h

Ei h

where the observed energy deposition

rate at altitude i is and the modeled one is Ei . For the functions that h compare the integrated energy deposited, we use the trapezoidal rule to define the difference between the cumulative energy deposited until altitude index i as

i

i 1

=

Ei 1 h

+

Ei h

E^i 1 h

Energy deposition rate functions tested:

E^i h

h 2

• Root-mean-square error (RMSE): min • Maximum error or runout: min (max(|ε |)) • Average error: min

n 2 i=1 i

n

1 n

n i=1

i

i

Cumulative energy deposition functions tested:

• Total energy error: min δ • RMSE of cumulative energy deposition: min • Runout of cumulative energy deposition: min (max(|δ |)) • Average error of cumulative energy deposition: min n

n 2 i=1 i

n

1 n

i

n i=1

i

Statistical association functions tested:

• Correlation: max corr ( , ) • Peaks and troughs matching: We also implemented an objective E h

E h

function that compared the number of peaks and troughs in addition to calculating the distances of peaks and troughs between the observed and modeled curves.

The objective function that consistently outperformed the others in capturing fragmentation height and energy was the runout of the energy deposition rate. Moreover, the objective function that exceeded the others in matching the curve's general shape was the RMSE of the energy deposition rate. When used together, these techniques ensure that the main trends of the curve are emulated by feature matching at the largest peaks and by minimizing the deviation through the entirety of the curve thereafter. The RMSE drives the solution in the initial generations but then becomes secondary as the GA evolves the fits. Typically, by the end of the evolution, the runout is an order of magnitude greater than the RMSE. GA solutions using this approach provided consistent curve matching across a range of synthetic curves and serve as the basis of the results in this study. The combined objective function used in this study is

Table 1 Ranges of the asteroid properties and FCM parameters that comprise the decision variables and bound the search space for the meteor modeling GA.

3

=

Ei h

2.3.1. Parameter search space In this study, the parameter values for each case within the initial population are selected at random from uniformly distributed userdesignated ranges. Table 1 displays the range of values for each of the decision variables. In addition to the initial asteroid properties of interest, we also treat the FCM fragmentation and ablation parameters (σ, α, number of fragments, and cloud fraction) as decision variables for the GA so that the optimizer can fit the best solution without the user overspecifying the modeling assumptions. However, these parameters are treated as modeling uncertainties and are not reported as part of the meteor inference in this study. The decision variables are continuous, positive, real numbers, with the exception of the number of fragments, which is an integer. The search space is large, especially with several variables spanning several orders of magnitude. Any of the genes can be set as a constant among all cases, which is especially helpful when some values are known. In this study, we have assumed published values for the initial entry angle and velocity. The FCM GA implementation and parameter search space in this study focus on larger initial meteoroid sizes, limited to diameters > 10 cm, to ensure that modeled bodies are in the continuum flow regime when they reach the initial point of the simulation at 100 km altitude. Smaller meteoroids require other modeling approaches, such as thermal fragmentation (e.g., Campbell-Brown and Koschny, 2004). Since our objective is to facilitate inference for asteroid threat assessment,

Parameter

i

273

Icarus 329 (2019) 270–281

A.M. Tárano, et al.

min 10

n 2 i=1 i

n

+ max(| i|)

Table 2 GA hyperparameters using steady-state population model.

(1)

where εi represents the vector containing the residuals of length n, representing the energy deposition rate at a given altitude i. This function is used to evaluate the fit of each FCM case in the population, with lower objective values indicating better fits (i.e., lower error compared to the observed curve). 2.3.3. Evolution process and genetic operators We tested a variety of GA selection mechanisms, including uniform, tournament, rank, and Roulette wheel, and several overlapping schemes, such as steady-state and generational GAs, as described in Wall (1996). We evaluated the performance of each implementation by the ability of the implementation to reproduce the synthetic curves and associated parameter values independent of random seed initiation, given a constant number of generations, population size, and genetic operators. The GA evolution approach that we found to be best suited to our meteor modeling cases uses a steady-state population model (Jones and Soule, 2006) with a Roulette wheel selection mechanism (Goldberg, 1989). In this approach, the population size remains fixed from generation to generation, and cases are selected for advancement and possible mating and mutation based on their fitness level (Wall, 1996). In the replacement process, the best performing individuals—limited by the prescribed population size—move on to the next generation, regardless of whether they are parent, child, or from the current generation. Controlled experiments were performed to determine the best-performing GA hyperparameters, including the population size, number of generations evaluated, mating probability, mutation probability, and the fraction of cases replaced in each iteration. A summary of the evolution process is presented in Fig. 2. We opted for the steady-state population model rather than a generational approach because of the increased convergence speed realized in our tests. However, with a steady-state GA we must ensure the population size is large enough as to avoid premature convergence to local minima. We tested the population sizes from 10,000 to 200,000, and determined that 80,000 individuals evolved for 800 generations provided similar objective scores within a reasonable computational time. The Roulette wheel selection mechanism stochastically selects pairs for reproduction based on their fitness level; the probability of being chosen for mating is proportional to the individual's fitness over that of the entire population (Goldberg, 1989). In our optimal design, 70% of the current population in each generation is selected for potential mating and mutating. These candidates are paired and each pair has a 75% probability of mating and undergoing genetic recombination. An even-and-odd crossover operator is used to mate the pairs selected for recombination. In this approach, two children are produced by alternating every other gene between the mother and father: one using genes with even indices from the mother and odd indices from the father, and vice versa for the second child. This crossover type allows for reproduction that explores diverse areas if the two parents have very dissimilar features. The remaining 25% of pairs not selected for mating

Parameter

Value

Number of genes Population size Number of generations Probability of crossover Probability of replacement Probability of mutation

7 80,000 800 0.75 0.70 0.72 (0–200 generations) 0.50 (201–400 generations) 0.15 (401–800 generations)

advance unchanged. Next, the genes of the advancing population (both mated and unchanged) are stochastically selected for potential mutation. Mutation introduces new genetic material into the pool through a random reselection of some of the genes, sampled from the same ranges given for the initial population (Table 1). In order to mitigate the potential for premature convergence to a local minimum, we utilized an adaptive mutation probability (Uyar et al., 2004) that decreases deterministically as the generations increase. A mutation probability of 0.72 is used for the first 200 generations, 0.5 for the next 200, and 0.15 for the remaining generations. The replacement mechanism is set so that only the best performing individuals—limited by the fixed population size—move on to the next generation, regardless of whether they are an original member of the original generation or a member of the newly generated population. Table 2 summarizes the parameters that will be used for the rest of the study. 3. Results The methodology described above will now be applied to three observed meteor cases: Chelyabinsk, Lost City, and Benešov. We chose these three events because they are well documented and present a variety of different fragmentation features, including a large primary flare, a smooth ablation-dominated flare, and multiple smaller flares. For each case, we ran the GA 15 times with different random seeds to provide an estimate of the variability of the approach due to randomness in the process. The resulting density, diameter, mass, and strength results are compared to published values to demonstrate the GA based inference. We will briefly summarize each event then present and discuss the results from the GA inference. 3.1. Chelyabinsk The superbolide that airburst over Chelyabinsk, Russia on February 15, 2013 is one of the largest and most well-documented meteor events available. The meteoroid was initially traveling 19.16 km/s at a flight path angle of 18.3° with respect to the horizontal (Popova et al., 2013). We use these trajectory parameters for this GA example. Many researchers have used this meteor for validation of simulations (e.g.,

Fig. 2. A typical genetic algorithm (GA) structure is presented from initiation to finding the best solution. The evolution process starts with an initial population. This population is then evaluated. While the stopping criterion isn't achieved, individuals are selected for potentially mating and mutating. The offspring and the parents are evaluated by the objective function so that only the best continue on to next iteration. 274

Icarus 329 (2019) 270–281

A.M. Tárano, et al.

Fig. 3. (a) A sample solution of the GA to the Chelyabinsk curve assuming the entry angle and velocity. This example shows that the two main peaks were captured. The parameters that led to this curve can be found in Table 3. The blue curve shows the best GA solution, and the gray dotted lines represent the range of the fittest individuals from all the GA runs. (b) The residual error from the best GA solution, plotted as the difference in energy deposition between the GA solution and the target meteor curve at each altitude increment (shown in linear scale).

Avramenko et al., 2014) and characterization of meteorite properties (e.g., Morlok et al., 2017). Using X-ray computed tomography, a meteorite-derived density of 3.3 g/cm3 was measured (Popova et al., 2013). The energy deposition curve used for this application was derived from the light curve synthesized from video observations (Brown et al., 2013). The best solution found by the GA and its error at each altitude are shown in Fig. 3. The GA was particularly capable of reproducing the altitudes and peak energy deposition rates of the two largest flares. The GA prioritizes the matching of these main peaks since they hold more weight in the objective function, from a magnitude consideration. The curvature of the main peak is well reproduced at the maximum energy deposition point of 80 kT/km at 30 km of altitude. Over the many GA runs we conducted, the GA generated solutions that match the main peak early in the evolution process. The GA then focused on fitting the secondary peak in the subsequent generations. Although, the secondary peak is narrower in the FCM solution than the observed curve, the altitude is within a few percent of the observed value. These results indicate that the objective function properly rewards solutions that match the most massive flares, which is critical for meaningful parameter inference. However, the GA was unable to capture the minor flare seen at approximately 41 km altitude because we assume a single monolithic body at the top of the atmosphere. Wheeler et al. (2017, 2018) found that Chelyabinsk's upper flare could not be reproduced with the embedded assumptions of uniform fragmentation and uniform structure found in the Wheeler et al. (2017) FCM version, which most closely resembles the version employed in this study. Future work will likely address such features. The GA, however, was capable of reproducing the most salient features of the curve, such as the two main peaks and the low energy deposition rates above and below the observed values. The FCM inputs that produced the solution presented in Fig. 3 are 16.35 m for initial diameter, 4.66 g/cm3 for bulk density, and 0.90 MPa for initial breakup strength. These parameters yield an initial mass of 1.07 E7 kg. Table 3 contains a summary of the parameters along with the range of values across the multiple random seeds. In our current best solution, the initial diameter of 16.4 m most

Table 3 Results of 15 GA runs for Chelyabinsk. These values provide the range of solutions for each parameter from the fittest individual from each of the GA runs. Parameter

Best case

Minimum

Maximum

Mean

Standard deviation

Objective score Initial mass (kg) Diameter (m) Density (g/cm3) Initial strength (MPa) Ablation coefficient (kg/J)

12.98 1.07E7 16.35 4.66 0.90 3.78E-9

12.98 1.04E7 15.99 3.21 0.90 1.91E-9

16.29 1.08E7 18.49 4.99 1.30 4.39E-9

13.95 1.06E7 17.23 4.00 0.99 3.08E-9

0.77 8.93E4 0.69 0.48 0.08 6.86E-10

closely matches the size estimates of Popova et al. (2013), which range from 15.2 to 24.4 m and were derived from a light curve using video observations. Diameter estimates below the mean value corresponded to bulk densities over 4 g/cm3. The results with higher diameters and similar objective values to the best solution match the 18–19 m estimates from government satellites by Brown et al. (2016). Other studies used similar ranges of 18–20 m for the initial size of the asteroid even within different modeling techniques (Avramenko et al., 2014; Borovička et al., 2015; Wheeler et al., 2017; Wheeler et al., 2018). These higher diameter estimates corresponded to results where the bulk density is within 5% less than the calculated meteorite density of 3.3 g/ cm3 (Popova et al., 2013). This example shows the uncertainty that continues to exist within mass trade-offs. Published initial mass estimates range from 1 E7 kg to 1.3 E7 kg (Avramenko et al., 2014; Borovička et al., 2015; Brown et al., 2016; Popova et al., 2013; Wheeler et al., 2017, 2018), with a potential factor of two in uncertainty due to uncertain energy estimates (Popova et al., 2013). Our estimate of 1.1 E7 kg falls well within expected range for mass. The initial strength of 0.90 MPa the GA solves is most similar to the 1.3–1.6 MPa estimate from Wheeler et al. (2017), with fragmentation starting at similar altitudes below 40 km. Wheeler et al. (2017) were not 275

Icarus 329 (2019) 270–281

A.M. Tárano, et al.

able to reproduce the upper flare in their best fits due to the monolithic structure. The Weibull size-strength scaling is unable to reproduce the variability of flare sizes seen in actual meteor breakup since it means that smaller pieces are always stronger and break lower than larger pieces. In reality, smaller pieces may break off before larger primary disruptions take place. Therefore, different rubble boulders with a range of independent strengths and sizes are required for larger fragments to persist at lower altitudes to better represent multiple bursts, as suggested in Wheeler et al. (2017). Wheeler et al. (2018) have addressed the upper flare by considering initial asteroid structure other than a pure monolith. Avramenko et al. (2014), Borovička et al., (2015), and Popova et al. (2013) have significantly lower values of initial strength 0.2–0.5 MPa because they begin their fragmentation sequence around 45 km, where the minor flare is missed by the GA. As a result, the difference in strength between our current result and published values is expected. More discussion on this is available in the Discussion section. The current inferred initial mass and diameter values agree with those from the literature for the Chelyabinsk event. The large, broad primary flare is well suited for modeling. To more extensively test the approach, events with fundamentally different behaviors are examined in the following subsections.

In order to match the more specific debris-shedding or ablation variations in future work, more variable fragmentation modeling parameters would be needed, such as the extended version of FCM by Wheeler et al. (2018). The FCM inputs that produced the solution presented in Fig. 4 are 0.39 m for initial diameter, 4.97 g/cm3 for bulk density, and 7.82 MPa for initial breakup strength, which is strong enough that it doesn't fragment. Table 4 contains a summary of the inferred initial values. The initial diameters inferred by the GA were consistently within 0.38–0.40 m, which is at most 16% smaller than the 0.45-m estimate published by Popova et al. (2011). That estimate was based on an initial mass estimate of 160 kg and a bulk density derived from a class average of 3.4 ± 0.18 g/cm3 from Britt and Consolmagno (2003, 2004). Our GA solved for much higher bulk densities ranging from 4.8 to 5 g/cm3 assuming a spherical volume. When comparing the inferred initial masses from different GA runs, the Lost City estimates were all within 6% of each other, ranging from 150 to 160 kg, with 153 kg as the point estimate derived from the best fit. This result is consistent with other published mass estimates, which range from 149 to 171 kg (Ceplecha, 1996; Ceplecha and ReVelle, 2005; Popova et al., 2011). The GA found solutions with high initial strengths and high ablation coefficients to produce non-breaking ablation-driven cases. According to Mathias et al. (2017), the ablation coefficient parameter is highly uncertain, ranging from 3.5E-10 to 7E-8 kg/J. The GA results for Lost City continuously reported values of around 5E-8 kg/J, which is within the high side of this uncertainty range. Other Lost City modeling efforts, such as Ceplecha (1996), Ceplecha et al. (1993), and Stulov et al. (1995) (as cited in Gritsevich, 2008), have also computed moderately high ablation coefficient values of around 1.0E-8 to 1.5E-8 kg/J (0.01–0.015 s2/km2), which are on the same order of magnitude as the GA estimates. As discussed in Ceplecha and ReVelle (2005), these higher ablation coefficient values represent “apparent” ablation rates, which include contributions from fragmentation processes that are not explicitly treated in single-body models. Ceplecha and ReVelle (2005) also computed a lower “intrinsic” ablation coefficient estimate of 8E9 kg/J (0.008 s2/km2) for Lost City using their FM model, which does explicitly include fragmentation processes. However, since the GA model's best matches for Lost City were all non-breaking, single-body cases, we expect the ablation coefficient results to be closer to the higher “apparent” values. Our current approach was able to reproduce the primary feature of the light curve, and reasonably infer initial parameters to well within uncertainties from the published results. However the smaller, highaltitude flares and the recognized fragmentation events near 40 km and 22 km altitudes by Ceplecha (1996) were missed. Capabilities for matching multiple flares that are not as strongly pronounced remain the subject of future work.

3.2. Lost City The Lost City meteor event took place in the United States on January 3, 1970. The meteor was recorded by the Smithsonian Prairie Network and the meteorites were chemically analyzed (Clarke et al., 1971; Nava et al., 1971). The entry velocity and angle, with respect to the horizontal, were calculated to be 14.2 km/s and 38∘, respectively (Ceplecha, 1996; McCrosky et al., 1971; ReVelle and Ceplecha, 2002). Popova et al. (2011) summarized these values with other orbital origin information and meteorite analysis data. The meteorite analysis led to a material density calculation of 3.7 g/cm3 (Clarke et al., 1971; McCrosky et al., 1971). We derived the energy deposition curve from Lost City's light curve using the same methodology that Wheeler et al. (2018) applied for Benešov. We used the light curve, deceleration profile, initial trajectory information, and initial mass estimate of 165 kg from Ceplecha and ReVelle (2005). First, we computed the power radiated by the meteor as a function of time by assuming that a zero-magnitude meteor emits 1530 W of bolometric power in the panchromatic passband (Ceplecha and ReVelle, 2005). Then we integrated the power radiated over the entire trajectory to estimate the total energy radiated by the meteor. The initial kinetic energy of the meteoroid was computed using initial mass and speed from Ceplecha and ReVelle (2005). A constant luminous efficiency of 1.79% was calculated by taking the ratio of total radiated energy to initial kinetic energy. The energy deposition per unit time was computed by dividing the radiated power by the luminous efficiency. Finally, the energy deposition per unit height was found by dividing the deposition per unit time by the vertical speed of the meteoroid. The best solution found by the GA and its error at each altitude are shown in Fig. 4. The GA solutions consistently produced results where the general trend and width of the curve are well matched. Given the current approach, the gradual, gentle curve is best approximated using non-fragmenting, ablation-driven cases, since it did not exhibit any distinct flares from catastrophic disruptions. As a result, the gradual curvature and slope of the profile are well reproduced, but more smoothly. However, the upper altitude oscillations were never fully replicated by any of the GA solutions due to the assumption of a single monolithic body, similarly to the upper flare missed in the Chelyabinsk light curve. Two fragmentation events reported by Ceplecha (1996) near 40 km and 22 km altitudes were not captured either because the objective function rewards high ablation coefficients that capture the curve's general shape rather than the less defined flares from this curve.

3.3. Benešov The Benešov superbolide occurred May 7th, 1991 over the Czech Republic where a rich spectrum was recorded (Borovička and Spurnỳ, 1996). Twenty years later after a reanalysis of the trajectory, an expedition successfully recovered some fragments of the fall that demonstrated that the meteoroid was heterogeneous with meteorites of types LL3.5, H5, and primitive achondrite (Spurnỳ et al., 2014). Unfortunately, no bulk density measurements are available but the mean meteorite bulk densities of LL and H type falls imply a range from 3.2 g/ cm3 to 3.4 g/cm3, respectively (Flynn et al., 2018; Macke et al., 2010). The trajectory calculations estimate that the meteoroid traveled 21.3 km/s with an entry angle of 81∘ as it entered the atmosphere (Borovička and Berezhnoy, 2016). We assume these initial entry trajectory parameters in the GA. The energy deposition curve that we used for Benešov is from Wheeler et al. (2018), and was based on the light curve and 4100-kg mass estimate of Ceplecha and ReVelle (2005). Fig. 5 shows the solution with the lowest objective value found by 276

Icarus 329 (2019) 270–281

A.M. Tárano, et al.

Fig. 4. (a) A sample solution of the GA to the Lost City curve assuming the entry angle and velocity. This example shows the GA's ability to emulate an energy deposition curve with gentle and gradual slope. The parameters that led to this curve can be found in Table 4. The gray dotted lines represent the range of the fittest individuals from all the GA runs. The range is indistinguishably different from the best result. (b) The residual error from the best GA solution, plotted as the difference in energy deposition between the the GA solution and the target meteor curve at each altitude increment (shown in linear scale).

for solutions that provided initial diameters between 1 and 1.3 m. The higher diameter values are closer to the diameter estimates of most modelers, who used 1.3 and 1.35 m (Brown et al. (2016); Wheeler et al., (2018), respectively). The bulk density from the best GA run is well within the 2–3.7 g/cm3 range used by other modelers (Borovička and Spurnỳ, 1996; Brown et al., 2016; Wheeler et al., 2018). The GA runs with higher diameters were paired with lower densities (2.7–3.3 g/cm3) than the best fit presented in Fig. 5. These solutions had objective values less than a percent different from that of the best case, exemplifying the typical mass trade-offs that lead to similar curves. However, the GA provides a set of potential solutions within a large parameter space. Within all the GA runs performed for this example, the initial mass estimates were within < 10% of 2706 kg. These values are all within the initial mass estimates of other studies, which have ranged from 75 kg to 4200 kg, with most estimates between around 2000–4000 kg. Those prior estimates are based on a range of modeling techniques, such as dynamically derived values and radiative hydrodynamic modeling (Borovička and Berezhnoy, 2016; Borovička et al., 1998; Borovička et al., 2015; Brown et al., 2016; Ceplecha and ReVelle, 2005; Spurnỳ et al., 2014; Wheeler et al., 2018). Compared to the 4100 kg estimate used to derive the Benešov energy deposition curve, our solutions underestimate the initial masses because the width of various peaks are not well reproduced and the top section of the curve is missed. This result suggests that the current GA design may underestimate the initial mass of the meteoroid if there are many peaks of significant width within the curve. The initial breakup fragmentation event of the GA solution is at a lower altitude than other published modeling fits, which typically started fragmentation well above 34 km (Borovička and Spurnỳ, 1996; Borovička et al., 1998; Borovička et al., 2015; Wheeler et al., 2018). Therefore, the initial breakup strengths reported from external sources are much smaller–some being 0.02 MPa, corresponding to breakup above 60 km of altitude (Wheeler et al., 2018). When comparing the strengths of subsequent breakup events associated with the main peaks rather than the initial breakup, however, the GA strength value agrees

Table 4 Results of 15 GA runs for Lost City. These values provide the range of solutions for each parameter from the fittest individual from each of the GA runs. Parameter

Best case

Minimum

Maximum

Mean

Standard deviation

Objective score Initial mass (kg) Diameter (m) Density (g/cm3) Ablation coefficient (kg/J)

3.32E-4 153 0.39 4.97 5.14E-8

3.32E-4 151 0.39 4.97 5.09E-8

3.33E-4 153 0.39 5.00 5.15E-8

3.32E-4 152 0.39 4.99 5.12E-8

3.80E-7 0.47 5.51E-4 8.14E-3 1.63E-10

the GA and its error at each altitude when assuming the initial trajectory for Benešov. This result demonstrates the GA's ability to capture the altitudes and approximate magnitudes for multiple energy deposition peaks. The locations and the values of the two main peaks are prioritized by the GA's objective function by minimizing error in both the RMSE and runout calculations. Furthermore, both peaks match the altitude of the observed flares with negligible difference. Although the maximum energy deposition rates at these peaks are overestimated by the GA solutions, the values are only about 2 times higher than the observed values, which is reasonable considering the uncertainty within the energy deposition measurements and models. These results were consistent among the multiple GA runs, and shown by the bounding gray lines in Fig. 5. The parameter inputs that produced the lowest objective value solution shown in Fig. 5 were 1.14 m for initial diameter, 3.44 g/cm3 for bulk density, and 3.76 MPa for initial breakup strength. These inferred initial parameters led to an estimated 2706 kg initial mass, assuming a spherical volume. The ranges of results from the multiple GA runs are summarized in Table 5. Overall, the GA runs robustly provided estimates of initial diameter and mass while the bulk density was not as well constrained. For the solution in Fig. 5, the initial diameter most closely resembles the Borovička et al. (2015) estimate of > 1 m. The GA consistently solved 277

Icarus 329 (2019) 270–281

A.M. Tárano, et al.

Fig. 5. (a) A sample solution of the GA to the Benešov curve assuming the entry angle and velocity. This example shows that the two main peaks were captured. The parameters that led to this curve can be found in Table 5. The gray dotted lines represent the range of the fittest individuals from all the GA runs. (b) The residual error from the best GA solution, plotted as the difference in energy deposition between the GA solution and the target meteor curve at each altitude increment (shown in linear scale).

search space through simplified modeling assumptions to decrease computation time and number of local minima. Although Borovivcka (2015) found no evidence that meteoroids 1–20 m in size are likely to be rubble piles, modelers such as Wheeler et al. (2018) or Borovička et al. (2013) have needed to employ heterogeneous structures or breakup mechanisms to capture the multiple disruptions and diverse flare features exhibited by many observed meteors. To enable these more detailed matches and inferences, the GA could be enhanced to include more complex fragmentation modeling, with additional decision variables. Future work will extend the GA approach to include the rubble pile and fragmented body modeling extensions from Wheeler et al. (2018). Using this more complex version of FCM, Wheeler et al. (2018) were able to obtain more detailed manual fits to the energy deposition curves for Chelyabinsk and Benešov than is possible using the more simplified fragmentation assumptions employed for our initial GA development. In those matches, Wheeler et al. (2018) assumed the published estimates for the asteroid properties, and used the more complex fragmentation approach to make inferences about the initial structure and breakup characteristics. In this development study, we have instead begun by assuming more simplified fragmentation parameters and inferring the main initial asteroid properties. The extended FCM version in Wheeler et al. (2018) allowed for a different ablation coefficient for clouds and fragments, uneven fragment mass-split distributions, a variable cloud dispersion coefficient, and a non-monolithic representation of the entering body. This approach allows fragment groups of various sizes and properties to begin disruptions at different altitudes, which enables fitting of a broader range of flare features that are not well represented by Weibull-based size-strength scaling or even fragment splits. For example, Chelyabinsk's upper hump was fit using a non-monolithic body, which required either a debris cloud or a release of a small, weak fragment with high cloud mass fraction. Similarly, the Benešov match in Wheeler et al. (2018) required many small fragments with different strength values to reproduce the gradual breakup along the upper curve and the sequence of three similar-sized sub-flares in the nose of the peak. Lost City was

Table 5 Results of 15 GA runs for Benešov. These values provide the range of solutions for each parameter from the fittest individual from each of the GA runs. Parameter

Best case

Minimum

Maximum

Mean

Standard deviation

Objective score Initial mass (kg) Diameter (m) Density (g/cm3) Initial strength (MPa) Ablation coefficient (kg/J)

0.058 2706 1.14 3.44 3.76

0.058 2616 1.03 2.79 3.53

0.059 2897 1.26 4.81 4.06

0.058 2764 1.13 3.70 3.84

3.18E-4 83.24 0.068 0.63 0.17

2.72E-8

2.29E-8

3.68E-8

3.05E-8

4.32E-9

quite well with those reported by Borovička and Spurnỳ (1996), Borovička et al. (2015), and Wheeler et al. (2018). For example, Borovička and Spurnỳ (1996) found best matches with breakup strengths of 4 MPa at 33 km and Borovička et al. (2015) of 2.5 MPa at 38 km. The initial strength solutions for the many GA runs were well constrained within 3.5–4.1 MPa. 4. Discussion The GA approach could be extended to be used with a range of other entry and fragmentation models—from simpler progressive fragmentation or pancake models to more complex models, such as the rubble pile extension of FCM presented in Wheeler et al. (2018)—depending on the given characterization objectives. Simpler models provide an efficient means of approximating basic properties or average characteristics, but lack the fidelity to match specific features that may offer more detailed inference. More complex models have the capability to match more detailed features, but also introduce more free parameters that increase the potential search space. Particular care must be taken to ensure that an increase in search space is paralleled with a GA design that proactively encourages sufficient diversity through the iterations. The FCM version used in this paper reduces the dimensions of the 278

Icarus 329 (2019) 270–281

A.M. Tárano, et al.

not modeled by Wheeler et al. (2018). However, we expect that using the extended FCM version may allow for the recorded fragmentation events to be captured using our approach. We predict that distinct ablation coefficients and different strength values assigned to different rubble groups could capture the widths of the flares that are not as strongly pronounced without requiring ablation to drive the energy deposition. For both Chelyabinsk and Benešov, Wheeler et al. (2018) note that reducing the cloud dispersion coefficient facilitates the match of the broader width of the flares. This paper's version of the model uses the baseline cloud dispersion coefficient of 3.5 from Hills and Goda (1998) along with identical fragment splits and Weibull-based strength scaling, and therefore produces narrower flares with a constrained cascade of magnitudes. For Benešov, the GA compensated for this by finding solutions with low cloud mass fractions and high strength scaling parameters to produce sharp peaks located at the appropriate altitudes. However, the GA was able to successfully find the best match possible for the given model simplifications, and was able to obtain reasonable estimates for the main initial properties despite the more approximate energy depositions matches. We obtained similar initial mass estimates as those in Wheeler et al. (2017, 2018) even though our diameter and bulk density estimates differed. The GA approach allows for a mapping of the correlated properties and parameters within the search space. However, it does not solve the inherent non-uniqueness issue that occurs in meteor modeling (e.g., Borovička et al., 1998; Popova et al., 2013) due to energy and mass trade-offs or numerous free parameters. Rather, the GA provides a quantifiable fitness measure to compare the different possible solutions and to determine the possible range of those modeling parameters that are not well constrained. For this initial development effort, we have circumvented the energy trade-off issue by assuming the initial entry angle and speed. Even with the synthetic curve tests used to determine the GA hyperparameters, energy and mass trade-offs produced different sets of parameter values but similar energy deposition curves, with only the bottom and top of curves differing notably. However, a similar mass trade-off remains between density and diameter. For both Chelyabinsk and Lost City, the GA bulk density estimates were consistently higher than the meteorite measurements while maintaining similar initial mass inferences. The fits tend to match the energy (and by extension the mass when velocity is fixed), while they are less sensitive to trade-offs between diameter and density that produce similar masses and energies. In test cases where we assumed lower fixed densities for Chelyabinsk and Lost City—using one-fourth of individuals and half of the generations from the Methodology section, which are values determined viable from synthetic curve experiments—the resulting objective values were slightly worse than the mean objective values presented in Tables 3 and 4. For Chelyabinsk, the lower density solutions produced secondary peaks that were higher and narrower than those of the low-density solutions, which led to higher objective error values because of the effect on the runout component. Modeling factors that contribute to this effect include the even fragment mass-split simplification and the assumed debris cloud spread rate, which produce narrower peaks due to multiple small, fast-spreading debris clouds released at identical altitudes. Given the nominal fixed cloud dispersion coefficient value of 3.5 adopted in the cloud dispersion formula for this study (Fig. 1), higher density values yield slower-spreading clouds and better objective value scores in the current GA implementation. Incorporation of improved spread rate models, variable dispersion coefficient values, or uneven fragment splits may help to resolve this high-density bias. For Lost City—which undergoes a more gradual disintegration that is best represented by non-breaking ablation-driven cases in our model—the GA favored the smaller, higher-density mass-equivalents because they pushed the nose of the curve down slightly, bringing it closer to the observed peak at the very lower edge of the flare and thus minimizing the runout error. Because of the constant ablation

coefficient assumed in the current mass-loss formula (shown in Fig. 1), smaller initial cross-sectional areas serve to reduce the ablation early on, delaying more of the mass-loss to lower altitudes with higher air densities. However, ablation rates are expected to vary with altitude, velocity, area, and other properties (Johnston et al., 2018). Incorporation of such variable ablation rates throughout the entry may improve inference capabilities for smaller more ablation-driven cases, improve matches at high altitudes, and facilitate capturing the subtle flares corresponding to fragmentation. Similarly, modeling light emission with variable luminous efficiency rates, rather than assuming a constant average value to translate light curves to energy deposition, may also be beneficial. In these ways, even when correlated physical parameters or modeling simplifications make unique inference difficult, the GA results can be useful for revealing the range of possible properties and/or identifying modeling simplifications or assumptions that could benefit from further refinement. 5. Conclusions and future work In this paper, we examined the capacity of a genetic-algorithmbased optimizer to automate meteoroid parameter inference using FCM. The GA parameters were optimized using an FCM-derived set of curves, then the resulting model was used to infer initial values of diameter, density, strength, and mass for three test cases: Chelyabinsk, Lost City, and Benešov. The approach is able to predict initial values, to well within published uncertainty, for a variety of events with different energy deposition characteristics. Moreover, once the GA hyperparameters are set, the inference can be performed automatically without extensive user-in-the-loop activity. This ability opens the potential for assessing a large number of events, which is critical to refining asteroid property distributions and uncertain modeling parameters for impact risk assessment. This methodology is also applicable to the wider asteroid and meteor community, who depend on modeling to provide inference for specific observed events. In the future, we will update the model to include the FCM rubble pile capabilities from Wheeler et al. (2018). In addition, we will adapt the objective function to directly compare light curves, eliminating the mapping to energy deposition profiles. Finally, with the updated approach, we will consider cases without a known entry trajectory to extend the ability to assess a much larger sample of bolide events. Acknowledgments This work was supported by the NASA Planetary Defense Coordination Office under contract NNA16BD60C. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. The software for this work used the GAlib genetic algorithm package, written by Matthew Wall at the Massachusetts Institute of Technology. The authors thank Jacqueline Machesky for her assistance in generating GA solutions. References Artemieva, N.A., Shuvalov, V.V., 2001. Motion of a fragmented meteoroid through the planetary atmosphere. Journal of Geophysical Research: Planets 1060 (E2), 3297–3309. Avramenko, Mikhail I., Glazyrin, Igor V., Ionov, Gennady V., Karpeev, Artem V., 2014. Simulation of the airwave caused by the Chelyabinsk superbolide. Journal of Geophysical Research: Atmospheres 1190 (12), 0 7035–7050. https://doi.org/10. 1002/2013JD021028. Baldwin, Barrett, Sheaffer, Yvonne, 1971. Ablation and breakup of large meteoroids during atmospheric entry. J. Geophys. Res. 760 (19), 0 4653–4668. Blum, Christian, Roli, Andrea, September 2003. Metaheuristics in combinatorial optimization: overview and conceptual comparison. ACM Comput. Surv. (0360-0300) 350 (3), 268–308. https://doi.org/10.1145/937503.937505. Borovička, Jiří, Berezhnoy, Alexey A., 2016. Radiation of molecules in Benešov bolide spectra. Icarus (0019-1035) 278 (0), 248–265.

279

Icarus 329 (2019) 270–281

A.M. Tárano, et al. Borovička, Jiří, Spurnỳ, Pavel, 1996. Radiation study of two very bright terrestrial bolides and an application to the comet s–l 9 collision with Jupiter. Icarus 1210 (2), 484–510. Borovička, Jiří, Popova, Olga P., Nemtchinov, Ivan V., Spurnỳ, Pavel, Ceplecha, Zdeněk, 1998. Bolides produced by impacts of large meteoroids into the earth's atmosphere: comparison of theory with observations. I. Benesov bolide dynamics and fragmentation. Astron. Astrophys. 334, 0 713–728. Borovička, JiřÍ, Tòth, Juraj, Igaz, Antal, Spurnỳ, Pavel, Kalenda, Pavel, Haloda, Jakub, Svoreň, Ján, KornoÅ¡, Leonard, Silber, Elizabeth, Brown, Peter, Marek, Husárik, 2013. The Košice meteorite fall: atmospheric trajectory, fragmentation, and orbit. Meteorit. Planet. Sci. 480 (10), 0 1757–1779. Borovička, Jiří, Spurnỳ, Pavel, Brown, Peter G., 2015. Small near-earth asteroids as a source of meteorites. In: Asteroids IV, pp. 257–280. Borovička, Jiří, Spurnỳ, Pavel, Grigore, Valentin I., Svoreň, Ján, 2017. The January 7, 2015, superbolide over Romania and structural diversity of meter-sized asteroids. Planetary and Space Science (0032-0633) 143, 0 147–158 (SI:Meteoroids 2016). Britt, Daniel T., Consolmagno, Guy J., 2003. Stony meteorite porosities and densities: a review of the data through 2001. Meteorit. Planet. Sci. 380 (8), 1161–1180. Britt, Daniel T., Consolmagno, Guy J., 2004. Meteorite porosities and densities: a review of trends in the data. In: Lunar and Planetary Science Conference. 35. Brown, Peter G., Revelle, Douglas O., Tagliaferri, Edward, Hildebrand, Alan R., 2002. An entry model for the Tagish lake fireball using seismic, satellite and infrasound records. Meteorit. Planet. Sci. 370 (5), 0 661–675. Brown, Peter G., Assink, Jelle D., Astiz, Luciana, Blaauw, Rhiannon, Boslough, Mark B., Borovička, Jiří, Brachet, N., Brown, D., Campbell-Brown, M., Ceranna, L., Cooke, William J., de Groot-Hedlin, C., Drob, D.P., Edwards, W., Evers, L.G., Garces, M., Gill, J., Hedlin, M., Kingery, A., Laske, G., Le Pichon, A., Mialle, P., Moser, D.E., Saffer, A., Silber, Elizabeth, Smets, P., Spalding, Richard E., Spurný, Pavel, Tagliaferri, Ed, Uren, D., Weryk, R.J., Whitaker, R., Krzeminski, Zbigniew, 2013. A 500-kiloton airburst over Chelyabinsk and an enhanced hazard from small impactors. Nature 503, 238 EP –, 11. Brown, Peter G., Wiegert, Paul, Clark, David L., Tagliaferri, Ed, 2016. Orbital and physical characteristics of meter-scale impactors from airburst observations. Icarus (00191035) 266, 96–111. Campbell-Brown, M.D., Koschny, D., 2004. Model of the ablation of faint meteors. Astronomy & Astrophysics 4180 (2), 751–758. Ceplecha, Zdeněk, 1996. Luminous efficiency based on photographic observations of the lost city fireball and implications for the influx of interplanetary bodies onto earth. Astron. Astrophys. 311, 0 329–332. Ceplecha, Zdeněk, ReVelle, Douglas O., 2005. Fragmentation model of meteoroid motion, mass loss, and radiation in the atmosphere. Meteorit. Planet. Sci. 400 (1), 0 35–54. https://doi.org/10.1111/j.1945-5100.2005.tb00363.x. Ceplecha, Zdeněk, Spurnỳ, Pavel, Borovička, Jiří, Keclíková, J., 1993. Atmospheric fragmentation of meteoriods. Astron. Astrophys. 279, 615–626. Ceplecha, Zdeněk, Borovička, Jiří, Graham Elford, W., ReVelle, Douglas O., Hawkes, Robert L., Porubčan, VladimÍr, Šimek, Miloš, Sep 1998. Meteor phenomena and bodies. Space Sci. Rev. (1572-9672) 840 (3), 327–471. Chyba, Christopher F., Thomas, Paul J., Zahnle, Kevin J., 1993. The 1908 Tunguska explosion: atmospheric disruption of a stony asteroid. Nature 361 0 40 EP –, 01. Clarke, Roy S., Jarosewich, Eugene, Nelen, Joseph, 1971. The Lost city, Oklahoma, meteorite: an introduction to its laboratory investigation and comparisons with pbram and ucera. J. Geophys. Res. 760 (17), 4135–4143. https://doi.org/10.1029/ JB076i017p04135. Cooke, William J., Moser, Danielle E., 2012. The status of the nasa all sky fireball network. In: Proceedings of the International Meteor Conference, 30th IMC, Sibiu, Romania, 2011, pp. 9–12. Demiriz, Ayhan, Bennett, Kristin, Embrechts, Mark J., 1999. Semi-supervised clustering using genetic algorithms. In: In Artificial Neural Networks in Engineering. ASME Press, pp. 809–814 ANNIE-99. Flynn, George J., Consolmagno, Guy J., Brown, Peter G., Macke, Robert J., 2018. Physical properties of the stone meteorites: implications for the properties of their parent bodies. Chemie der Erde - Geochemistry (0009-2819) 78 (3), 269–298. Gen, Mitsuo, Cheng, Runwei, 1997. Genetic Algorithms and Engineering Design. A Wiley Interscience Publication Wiley9780471127413. Gi, Nayeob, 2017. Using Bolide Airwaves to Estimate Meteoroid Source Characteristics and Window Damage Potential. Master's thesis. The University of Western Ontario, Canada. Goldberg, David E., 1989. Genetic Algorithms in Search, Optimization and Machine Learning, 1st edition. Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA0201157675. Gritsevich, M.I., Oct 2008. The pribram, lost city, innisfree, and neuschwanstein falls: an analysis of the atmospheric trajectories. Sol. Syst. Res. (1608-3423) 420 (5), 372–390. Hills, Jack G., Goda, M. Patrick, 1993. The fragmentation of small asteroids in the atmosphere. Astron. J. 105 (0), 1114–1144. Hills, Jack G., Goda, M. Patrick, 1998. Damage from the impacts of small asteroids. Planet. Space Sci. 460 (2), 219–229 ISSN 0032–0633. International Workshop Tunguska 96. Holland, John H., 1984. Genetic Algorithms and Adaptation, Pages 317–333. Springer US, Boston, MA978-1-4684-8941-5 ISBN. Huuk, Thiemo C., Hahn, Tobias, Osberghaus, Anna, Hubbuch, Jürgen, 2014. Model-based integrated optimization and evaluation of a multi-step ion exchange chromatography. Sep. Purif. Technol. 136, 207–222. ISSN 1383-5866. https://www. sciencedirect.com/science/article/pii/S1383586614005681. Johnston, Christopher O., Stern, Eric C., Wheeler, Lorien F., 2018. Radiative heating of large meteoroids during atmospheric entry. Icarus 309, 25–44 ISSN 0019-1035. Jones, Josh, Soule, Terry, 2006. Comparing genetic robustness in generational vs. steady

state evolutionary algorithms. In: Proceedings of the 8th Annual Conference on Genetic and Evolutionary Computation. ACM, New York, NY, USA, 1-59593-186-4, pp. 143–150. https://doi.org/10.1145/1143997.1144024. GECCO ’06. http:// doi.acm.org/10.1145/1143997.1144024. Levin, Boris I., 1956. The Physical Theory of Meteors, and Meteoric Matter in the Solar System. American Meteorological Society. Macke, Robert J., Consolmagno, Guy J., Britt, Daniel T., Hutson, Melinda L., 2010. Enstatite chondrite density, magnetic susceptibility, and porosity. Meteorit. Planet. Sci. 450 (9), 1513–1526. https://doi.org/10.1111/j.1945-5100.2010.01129.x. Mathias, Donovan L., Wheeler, Lorien F., Dotson, Jessie L., 2017. A probabilistic asteroid impact risk model: assessment of sub-300m impacts. Icarus (0019-1035) 289, 106–119. https://doi.org/10.1016/j.icarus.2017.02.009. McCrosky, Robert E., Posen, A., Schwartz, G., Shao, C.-Y., 1971. Lost city meteorite–its recovery and a comparison with other fireballs. J. Geophys. Res. 760 (17), 0 4090–4108. https://doi.org/10.1029/JB076i017p04090. Melosh, H.J., 1981. Atmospheric breakup of terrestrial impactors. In: Multi-Ring Basins: Formation and Evolution, pp. 29–35. Message Passing Interface Forum, 2015. MPI: A Message-Passing Interface Standard: Version 3.1. High-Performance Computing Center. https://books.google. com/books?id=uv1ajwEACAAJ. Mitchell, Melanie, 1996. An Introduction to Genetic Algorithms. MIT Press, Cambridge, MA, USA0-262-13316-4. Morlok, Andreas, Bischoff, Addi, Patzek, Markus, Sohn, Martin, Hiesinger, Harald, 2017. Chelyabinsk — a rock with many different (stony) faces: an infrared study. Icarus (0019-1035) 284, 0 431–442. Nava, David F., Walter, Louis S., Doan, Arthur S., 1971. Chemistry and mineralogy of the lost city meteorite. J. Geophys. Res. 760 (17), 0 4067–4071. https://doi.org/10. 1029/JB076i017p04067. Nemtchinov, Ivan V., Jacobs, C., Tagliaferri, E., 1997. Analysis of satellite observations of large meteoroid impacts. Ann. N. Y. Acad. Sci. 8220 (1), 0 303–317. Popova, Olga P., Jiří, Borovika, Hartmann, William K., Spurnỳ, Pavel, Gnos, Edwin, Nemtchinov, Ivan V., Josep, M.TrigoRodriguez., 2011. Very low strengths of interplanetary meteoroids and small asteroids. Meteorit. Planet. Sci. 460 (10), 0 1525–1550. https://doi.org/10.1111/j.1945-5100.2011.01247.x. Popova, Olga P., Jenniskens, Peter, Emel'yanenko, Vacheslav, Kartashova, Anna, Biryukov, Eugeny, Khaibrakhmanov, Sergey, Shuvalov, Valery, Rybnov, Yurij, Dudorov, Alexandr, Grokhovsky, Victor I., Badyukov, Dmitry D., Yin, Qing-Zhu, Gural, Peter S., Albers, Jim, Granvik, Mikael, Evers, Läslo G., Kuiper, Jacob, Kharlamov, Vladimir, Solovyov, Andrey, Rusakov, Yuri S., Korotkiy, Stanislav, Serdyuk, Ilya, Korochantsev, Alexander V., Larionov, Michail Yu., Glazachev, Dmitry, Mayer, Alexander E., Gisler, Galen, Gladkovsky, Sergei V., Wimpenny, Josh, Sanborn, Matthew E., Yamakawa, Akane, Verosub, Kenneth L., Rowland, Douglas J., Roeske, Sarah, Botto, Nicholas W., Friedrich, Jon M., Zolensky, Michael E., Le, Loan, Ross, Daniel, Ziegler, Karen, Nakamura, Tomoki, Ahn, Insu, Lee, Jong Ik, Zhou, Qin, Li, Xian-Hua, Li, Qiu-Li, Liu, Yu, Tang, Guo-Qiang, Hiroi, Takahiro, Sears, Derek, Weinstein, Ilya A., Vokhmintsev, Alexander S., Ishchenko, Alexei V., SchmittKopplin, Phillipe, Hertkorn, Norbert, Nagao, Keisuke, Haba, Makiko K., Komatsu, Mutsumi, Mikouchi, Takashi, 2013. Chelyabinsk airburst, damage assessment, meteorite recovery, and characterization. Science (0036-8075) 3420 (6162), 1069–1073. https://doi.org/10.1126/science.1242642. Pravec, Petr, Harris, Alan W., Kunirk, Peter, Gald, Adrin, Hornoch, Kamil, 2012. Absolute magnitudes of asteroids and a revision of asteroid albedo estimates from wise thermal observations. Icarus 2210 (1), 365–387 (ISSN 0019-1035). Register, Paul J., Mathias, Donovan L., Wheeler, Lorien F., 2017. Asteroid fragmentation approaches for modeling atmospheric energy deposition. Icarus (0019-1035) 284, 0 157–166. https://doi.org/10.1016/j.icarus.2016.11.020. ReVelle, Douglas O., 2007. NEO fireball diversity: energetics-based entry modeling and analysis techniques. Proceedings of the International Astronomical Union 20 (S236), 95–106. ReVelle, Douglas O., Ceplecha, Zdeněk, 2002. Semi-empirical fragmentation model of meteoroid motion and radiation during atmospheric penetration. In: Asteroids, Comets, and Meteors: ACM 2002. 500. pp. 285–288. Rumpf, Clemens M., Lewis, Hugh G., Atkinson, Peter M., 2017. Asteroid impact effects and their immediate hazards for human populations. Geophys. Res. Lett. 440 (8), 3433–3440. Seppelt, Ralf, Voinov, Alexey, 2002. Optimization methodology for land use patterns using spatially explicit landscape models. Ecol. Model. (0304-3800) 1510 (2), 0 125–142. Silber, Elizabeth A., Boslough, Mark, Hocking, Wayne K., Gritsevich, Maria, Whitaker, Rodney W., 2018. Physics of meteor generated shock waves in the earth's atmosphere — a review. Adv. Space Res. (0273-1177) 620 (3), 489–532. https://doi.org/10. 1016/j.asr.2018.05.010. http://www.sciencedirect.com/science/ article/pii/S027311771830406X. Spurnỳ, Pavel, Haloda, Jakub, Borovička, Jiří, Shrbenỳ, Lukáš, Halodová, Patricie, 2014. Reanalysis of the Benešov bolide and recovery of polymict breccia meteorites–old mystery solved after 20 years. Astronomy & Astrophysics 570 (0), A39. Stulov, V.P., Mirskii, V.N., Vislyi, A.I., 1995. Aerodynamics of Bolides. Science. Fizmatlit, Moscow. Stützle, T.G., 1999. Local Search Algorithms for Combinatorial Problems: Analysis, Improvements, and New Applications. Dissertationen zur künstlichen Intelligenz. Infix9783896012203. Tagliaferri, Edward, Spalding, Richard, Jacobs, Cliff, Worden, Simon P., Erlich, Adam, 1994. Detection of meteoroid impacts by optical sensors in earth orbit. In: Hazards Due to Comets and Asteroids. 24. pp. 199. Tsvetkov, V.I., Skripnik, A. Ia, 1991. Atmospheric fragmentation of meteorites from the standpoint of mechanical strength. Astronomicheskii Vestnik 25, 364–371.

280

Icarus 329 (2019) 270–281

A.M. Tárano, et al. Uyar, A.S., Eryigit, G., Sariel, Sanem, 2004. An adaptive mutation scheme in genetic algorithms for fastening the convergence to the optimum. In: Proceedings of 3rd APIS: Asian Pacific International Symposium on Information Technologies. Wall, Matthew, 1996. Galib: A C++ Library of Genetic Algorithm Components. 87. Mechanical Engineering Department, Massachusetts Institute of Technology, pp. 54. Wheeler, Lorien F., Register, Paul J., Mathias, Donovan L., 2017. A fragment-cloud model for asteroid breakup and atmospheric energy deposition. Icarus (0019-1035) 295,

149–169. https://doi.org/10.1016/j.icarus.2017.02.011. Wheeler, Lorien F., Mathias, Donovan L., Stokan, Edward, Brown, Peter G., 2018. Atmospheric energy deposition modeling and inference for varied meteoroid structures. Icarus (0019-1035) 315, 79–91. https://doi.org/10.1016/j.icarus.2018.06. 014. Zahnle, Kevin J., 1992. Airburst origin of dark shadows on venus. Journal of Geophysical Research: Planets 970 (E6), 10243–10255.

281