Performance of D-criteria in isolating meteor showers from the sporadic background in an optical data set

Separating meteor showers from the sporadic meteor background is critical for the study of both showers and the sporadic complex. The linkage of meteors to meteor showers, to parent bodies, and to other meteors is done using measures of orbital similarity. These measures often take the form of so-called D-parameters and are generally paired with some cutoff value within which two orbits are considered related. The appropriate cutoff value can depend on the size of the data-set (Southworth&Hawkins 1963), the sporadic contribution within the observed size range (Jopek 1995), or the inclination of the shower (Galligan 2001). If the goal is to minimize sporadic contamination of the extracted shower, the cutoff value should also reflect the strength of the shower compared to the local sporadic background. In this paper, we present a method for determining, on a per-shower basis, the orbital similarity cutoff value that corresponds to a chosen acceptable false-positive rate. This method also assists us in distinguishing which showers are significant within a set of data. We apply these methods to optical meteor observations from the NASA All-Sky and Southern Ontario Meteor Networks.


INTRODUCTION
Quantifying the similarity of meteoroid orbits to each other or to potential parent bodies is a critical first step in many meteoroid dynamics studies. Establishing a degree of orbital similarity between meteors is necessary in order to identify a new meteor shower, for instance. Orbital similarity is also presumed to exist between showers and their parent bodies. Finally, even the sporadic meteor complex largely consists of "sources" of orbitally similar meteors. Related bodies may be dissimilar due to meteoroid stream evolution, or they may simply appear dissimilar due to measurement uncertainties. Conversely, two meteors may resemble each other by chance. It is therefore important to establish the false-positive rate when using orbital similarity criteria to isolate showers or identify parent body candidates.
The degree to which two meteoroid orbits are similar is often evaluated using some formulation of a similarity parameter. Southworth & Hawkins (1963) formulated the first D-parameter, which is computed from two sets of orbital elements (see Section A1). D S H and its derivatives are larger for less similar orbits and are therefore sometimes referred to as "dissimilarity parameters" (e.g., Galligan 2001). Southworth & Hawkins (1963) designed their initial parameter to E-mail: althea.moorhead@nasa.gov be simple to compute and acknowledged that it is possible to formulate valid alternatives.
Alternatives were indeed formulated: one such variation is the Drummond D-parameter (D D , Section A2). The Drummond orbital similarity parameter (Drummond 1981) resembles D S H in its general form, but the four terms are scaled to have comparable weight and the last two terms use angular distances rather than chords. Jopek (1993) conducted a comparison of these two criteria and concluded that D S H is overly dependent on perihelion distance while D D is overly dependent on eccentricity. They proposed the use of a hybrid D parameter that inserts the pericenter term of D D into D S H . More recently, Valsecchi et al. (1999) introduced a new D-parameter (D N , Section A3) that is computed from observed meteor quantities: radiant, velocity, and time. D N is constructed such that two of its four components are nearinvariant under secular cycling of ω.
In order to use an orbital similarity parameter, some cutoff value must be adopted within which orbits can be considered related. Southworth & Hawkins (1963) modeled the false association rate between meteors by randomizing their nodes; they extracted meteor streams by linking meteors for which D S H ≤ 0.2. This value was specific to their dataset of 360 meteors; the critical value for single-linkage in larger or smaller data sets could be extrapolated from their value as 0.2 4 √ 360/N. Lindblad (1971) lowered this slightly to 0.8/ 4 Jopek (1995) later noted that the appropriate single-linkage cutoff is also a function of the relative percentage of sporadic meteors; Jopek & Froeschle (1997) computed linkage cutoffs by generating fake meteors using the overall meteor orbital element distributions. In each case, the single-linkage cutoff is a function of data-set size because it identifies meteors that are more closely related than average. Associating meteors with a known shower or parent body is more straightforward than performing cluster analyses. Generally, all shower meteors lying within a cutoff value are extracted and smaller datasets simply produce fewer extracted meteors. Drummond (1981) compared meteor showers with cometary parent bodies and found that good pairings generally fell within D S H ≤ 0.25 and D D ≤ 0.105. Similarly, Drummond (1991) suggested requiring D S H 0.2 to 0.25 and D D 0.1 to 0.125. Galligan (2001) tested the performance of D-criteria in recovering streams from a large radar data set, quoting cutoff values to recover 70% and 90% of a stream as a function of inclination. No matter how tight a cutoff is applied, however, there is some non-zero chance of including sporadic meteors.
Galligan modeled the level of sporadic intrusion in meteor stream recovery using two of the above orbital similarity parameters, D S H and D N . Figures 2 and 3 of Galligan (2001) display CDFs of sporadic AMOR meteors associated with stream searches as a function of D-parameter cutoff and inclination. We follow a similar approach; however, we present the sporadic intrusion for three D-parameters and for each individual shower considered rather than for inclination bands. We find it simpler to model the sporadic intrusion on a per-shower basis and compute the cutoff value for each shower at which the total contamination by sporadic meteors reaches 10%.
We apply our methods to data from the NASA All Sky Fireball Network (Cooke & Moser 2012) and Southern Ontario Meteor Network (Weryk et al. 2008). Combined, these systems have 25 all-sky meteor video cameras grouped in clusters in the U.S and Canada. Within a cluster, cameras are separated by 80-150 km; this configuration results in overlapping fields-of-view and permits meteor trajectory triangulation. Observations are automatically processed by the University of Western Ontario's ASGARD software (Weryk et al. 2008), which detects meteors and calculates trajectories. At the time of writing, these networks have measured trajectories for 27,113 meteors over the past 7 years.

METHODS
We investigate the performance of D S H , D D , and D N by comparing the shower association rate with a modeled falsepositive association rate (Section 2.1). We take the measurement uncertainties into account (Section 2.2), using them to assess the probability that a meteor belongs to a given shower. We investigate showers in order of decreasing significance in our data set, removing probable shower members with each iteration (Section 2.3).

False positive modeling
We compute the sporadic intrusion, or false positive rate for shower association, by analyzing the distribution of our orbital similarity parameters relative to shower "analogs" of our own construction. These analogs are positioned similarly in Sun-centered ecliptic coordinates, but are offset in time from the shower of interest by 60 • in solar longitude. These analogs are designed such that: [1.] they are offset from the original shower in time, [2.] the region of parameter space encapsulated by an orbital similarity cutoff is similar for these analogs as it is for the original shower orbit, and [3.] both shower and analogs have the same position relative to the sporadic sources in Sun-centered ecliptic coordinates. The only variations not accounted for by this method are seasonal variations in the sporadic sources and contributions from nearly showers. We attempt to minimize the influence of nearby showers by working our way from the strongest to the weakest shower, removing each along the way (see Section 2.3).
The exact construction of our shower analogs is as follows. We compute the shower's radiant in geocentric, Sun-centered ecliptic coordinates: λ g − λ and β g . We construct shower analogs by varying λ from λ ,shower + 60 • to λ ,shower + 300 • in increments of 10 • . The value of λ g varies at the same rate as λ so that λ g − λ is kept constant; β g and v g are also unchanged. The Sun-centered ecliptic radiant and geocentric velocity define the meteor stream's velocity relative to the Earth at the time of the shower. We compute the Earth's position and velocity at the given solar longitude in the year 2015 using the DE423 ephemeris. This allows us to complete the meteoroid stream's state vector at peak activity, which we then convert into orbital elements. This process yields the full set of shower parameters needed to compute D S H , D D , and D N .
For a perfectly circular Earth orbit, the above process would result in shower analogs that have the same orbital elements as the initial shower, aside from Ω, which would vary in the same manner as λ . We therefore model shower association false positives in an approach similar to that of Southworth & Hawkins (1963), who chose to randomize the node for their meteor data. We note, however, that due to the Earth's orbital eccentricity and the perturbations of the Moon, we produce analogs with values of q, e, i, and ω that are slightly different from the initial shower orbit.
The D-parameters for the shower analogs are combined into individual PDFs as well as a single CDF, the latter of which produces a smoothly increasing function that represents the number of expected false positives within a given cutoff value. Subtracting this predicted false positive rate from the shower D-parameter distribution yields an estimate of the true number of shower members within a given D value.

Uncertainty handling
Each of our meteor trajectory parameters has an estimated associated uncertainty, which can vary substantially between meteors. To limit the degree to which poorly-constrained orbits affect our results, we take the following measures. First, we apply some basic quality cuts. We require that the camera-meteor-camera angle Q * ≥ 15 • ; this geometry requirement favors better-determined trajectories. We also require that v g ≤ 75 km/s and that the associated uncertainty ∆v g ≤ (0.1 · v g + 1) km/s ( Figure 1). This dependence of ∆v g on v g mimics the observed correlation between the two quantities and is intended to avoid biasing the results against faster meteors.
After making these quality cuts, we create "clones" of our meteors by assuming Gaussian uncertainties in radiant and velocity (the uncertainty in solar longitude is effectively zero in comparison). We then compute the orbital elements corresponding to each ecliptic radiant and velocity choice. These clones are compared to both showers and shower analogs and each clone is individually evaluated for shower membership. We can then generate a non-binary shower membership probability for each meteor using the sum of its clones. This Monte-Carlo uncertainty handling produces smoother D-parameter distributions and allows us to compute shower membership probabilities that take measurement precision into account. This step can be omitted in order to analyze the shower association false-positive rate for meteor data that lacks estimated uncertainties.

Shower extraction
We developed a rough list of showers through the use of an orbital element heat map ( Figure 2). Showers were selected that correspond to local maxima in this map of ascending node and inclination. We also analyze several showers that are not visible in this heat map for the sake of demonstrating non-detections. For each shower, we perform our initial D calculations using previously established shower parameters (Table 1). We compute D-parameters for every meteor clone in our data set relative to our nominal shower orbit and radiant. If the number of low-D meteors exceeds the expected false positive rate, we use the ratio of the false positive D distribution to the shower D distribution to compute the probability that a meteor clone is not a member of the shower as a function of D. We then remove shower members using this distribution, repeating the process for the next strongest shower.
As we extract showers, we re-compute the shower's properties using its members. In each case, we compute the mean shower orbit as follows: where brackets indicate the arithmetic mean, the appropriate sign is taken for each arctangent computation, and the mean anti-velocity unit vector is − u = cos (λ − λ ) cos β x+ sin (λ − λ ) cos β ŷ + sin β ẑ. These values are converted to orbital elements, again using the DE423 ephemeris for the Earth to complete the state vector.

RESULTS
This study accomplishes two tasks. First, we compare the distribution of D-parameters about each shower with that about its analogs to determine whether a shower is detected. Second, we use the cumulative distribution of D-parameters to determine the maximum value of D within which the sporadic intrusion is no more than 10%.

Shower detection
Figure 3 compares D-parameter histograms for our data compared to the four most significant showers in our data set (black lines) and the false-positive rate computed from their analogs (gray region). For each shower, we present histograms for each of our three D-parameter choices. Plot ranges for each shower are synchronized, and each plot domain covers the same fraction of D i,max . We have also marked the proximity of other major showers in D-parameter space. The September -Perseids (SPEs), Orionids, η-Aquariids, and Leonids all lie within D D = 0.6 compared to the Perseids (PER). The influence of these showers is also reflected in the false positive rate computed using our Perseid analogs. Note from the middle-top panel of Figure 3 that the combined contribution of these showers lies within the predicted false positive rate.
Separating the Orionids from the η-Aquariids is more difficult. Neither D S H nor D D are able to distinguish between these two showers; the label marking the η-Aquariids in Figure 3 lies on top of the peak of the Orionids. Two peaks are visible in the D N distribution, but the η-Aquariids lie well above the predicted false positive distribution. This similarity is a result of the two showers being part of the same meteoroid stream -they should be grouped together by any measure of orbital similarity. To study these showers separately, it is necessary to slice the stream in two by, for instance, solar longitude.
The Northern and Southern Taurids (NTA and STA) are also challenging to separate using measures of orbital similarity; the only D-parameter capable of separating these branches is D S H . Unlike the Orionids and η-Aquariids, these two showers lie close to one another in both radiant and time of year. Our solution is to separate them by whether their radiant lies just above or just below the ecliptic. Meteor clone (N = 100) density map in longitude of ascending node (Ω) and inclination (i). Vertical striping occurs because Ω is generally measured more precisely than i. Blue areas represent low meteor density while red represents high density. The showers investigated in this paper are circled and labeled with the shower's three-letter code. Some are clearly visible, such as the Perseids (PER) and Geminids (GEM), while others appear minimal or absent, such as the epsilon Geminids (EGE) and Northern delta Aquariids (NDA).
plots that minor showers are overwhelmed by sporadics at lower D-parameter thresholds. Additionally, Figure 5 depicts two showers that we consider non-detections: the Northern δ-Aquariids (NDA) and the -Geminids. In both cases, the D N distribution is hardly distinguishable from our computed false positive rate.
We have also included two minor shower outbursts in this analysis: the May Camelopardalids (CAM) and the December Phoenicids (PHO). Both of these showers produced a tight cluster of five meteors in our data set that were automatically extracted from our data using the algorithm of Burt et al. (2014) and which corresponded to predicted outbursts. Figure 5 also indicates that each cluster of 4-5 meteors is anomalously similar compared to the modeled false positive rate. Figure 6 compares D-parameter cumulative distributions for three selected showers. The distribution about the nominal showers is depicted in black, while the distribution about the shower analogs, combined, appears as a solid gray line. By comparing these distributions, we can compute the cutoff value for each D-parameter at which the sporadic intrusion (shaded region) reaches 10% of the shower strength. Figure 6 also marks Galligan's 70% recovery criterion for each shower and D-parameter, where defined (Galligan 2001). For many of our major showers, including the Perseids, this cutoff is significantly smaller than our recommended cutoff. However, for the Leonids, Galligan's 70% recovery cutoff for D D is larger than our recommended cut-off. In this case, the blind use of Galligan's recommendation could result in including a large number of false positives.

Shower extraction
For smaller showers, such as the April Lyrids, Galligan's 70% cutoff values for each D-parameter produce very different numbers of shower members. These cutoffs would also produce significantly different false positive inclusion rates. This shower is a good illustration of the value of individually characterizing the false positive rate for a given shower and D-parameter.
In general, we agree with Galligan (2001) that D N parameter has the best overall performance. For most of our largest showers, D N retrieves the most shower members with 90% confidence (see Figure 7). For smaller showers, the three parameters vary in performance.

Shower orbits
We have analyzed the performance of each D-parameter relative to either established shower parameters, or, in a few cases, relative to shower orbits computed by the authors prior to beginning this study. By using previously measured orbits, we hope to avoid biasing our results in favor of one D-parameter over another. Now that we have completed our comparison of D-parameter performance, we compute average shower orbits using the method described in 2.3. Table 1 reports both the original reference orbits and our measured orbits. We report our computed shower orbits only for cases where the use of our computed orbit improves the shower removal (i.e., increases N 90 for D N ). In some cases, such as the Geminids, we were unable to improve on established orbits.  (Table 1). The shaded gray area encompasses the distribution of D parameters for each of our 25 "shower analogs" -the median appears as a solid gray line. The solid black line represents the difference between the raw D distribution for the shower and the median D distribution for the shower analogs. Cases where the D-parameter computed between the reference shower and another major shower falls within the depicted range are marked with small shower labels (e.g., "ETA"). In each plot, the dotted black line depicts the raw CDF of each D-parameter computed relative to the shower orbit (Table 1). The solid gray line marks the average CDF for the 25 "shower analogs." The solid black line represents the difference between the raw D cumulative distribution for the shower and the false positive D cumulative distribution for the shower analogs. The dashed blue vertical line marks Galligan's 70% shower retrieval cutoff for each D-parameter, where applicable, and shaded region marks the interval in which the sporadic intrusion is less than 10%.

CONCLUSIONS
We present a method for characterizing the expected false positive rate for shower extraction using orbital similarity criteria, or D-parameters. This method involves the construction of shower "analogs" that have the same Suncentered ecliptic radiant and geocentric velocity as the nominal shower, but vary in solar longitude. Shower detection occurs when the D-parameter distribution about the shower exceeds that about the shower analogs.
The construction of a false-positive distribution also permits the computation of the probability of shower membership as a function of D. It also assists the user in selecting an orbital similarity cutoff that limits the false positive rate to within a tolerable, user-determined percentage. We find that the traditional D-parameter cutoffs may fall short of or exceed what's appropriate, resulting in either needlessly throwing away shower members or including too many false positives. We recommend instead choosing a cutoff based on the modeled sporadic intrusion, especially for small show-ers which may be tricky to isolate from the sporadic background.
Like Galligan (2001), we find that the D N -criterion of Valsecchi et al. (1999) does the best overall job of isolating shower members from the sporadic background. For most of our strongest showers, D N was able to extract the largest number of meteors while limiting the sporadic intrusion to ≤ 10%.  Table 1. Reference orbits used to analyze the false-positive rate for D-parameter-based shower association. In a few cases, we generated our own reference orbit using meteors near peaks in an orbital element heat map ( Figure 2) or those selected automatically by a cluster detection algorithm (Burt et al., 2014). For each shower, we attempted to recalculate the average orbit using the method described in Section 2.3. We report these results only when extraction with the average orbit produced a higher yield than extraction with the original reference orbit.

A1 The Southworth & Hawkins D S H parameter
The first-established D-criterion that is still in wide use is that of Southworth & Hawkins (1963), which we will denote D S H . This parameter computes the degree of dissimilarity between two orbits as D 2 S H = (q 1 − q 2 ) 2 + (e 1 − e 2 ) 2 + 2 sin where I is the angle between the two orbital planes and Π is the angle between the two eccentricity vectors. These angles are computed as follows: In these equations, q j is the perihelion distance of meteoroid j, e j is its eccentricity, i j is its inclination, ω j is its argument of pericenter, and Ω j is its longitude of ascending node. The plus sign is used when Ω 1 and Ω 2 differ by less than 180 • , and the minus sign when the difference in node is greater than 180 • . Note that for non-hyperbolic, Earth-intersecting orbits, the first two terms of Equation A1 have a maximum value of unity, while the latter two terms have a maximum value of 4. Thus, D S H ≤ √ 10.
A2 The Drummond D D parameter (Drummond 1981) revised this criterion by substituting relative for absolute differences between orbital elements and by using angles in the place of chords (Galligan 2001). This revised parameter is Here, Π has been replaced by θ, which is defined as θ = arccos[sin β 1 sin β 2 + cos β 1 cos β 2 cos(λ 2 − λ 1 )] .
The use of fractional differences in the first two terms ensures that each of these terms cannot exceed unity. For nonhyperbolic orbits, D D ≤ 2.

A3 The Valsecchi et al. D N parameter
Valsecchi et al. (1999) proposed a new orbital similarity parameter, D N , that is more closely linked to observable quantities: geocentric right ascension and declination (α g and δ g ), geocentric velocity (v g ), and solar longitude (λ ⊕ ). These observables are converted into variables u, θ, and φ, where u is the ratio of v g to the Earth's velocity (v ⊕ ), θ describes the angle between v g and v ⊕ , and φ describes the component of the meteoroid's geocentric velocity that is perpendicular to v ⊕ .
These quantities can be calculated from observed quantities as follows.