HelioLinC NEO Study

I mentioned in my first HelioLinC study that I had found some new data to explore. Catalina Sky Survey has started pushing their nightly data to the Planetary Data System. CSS' mission is to find Near Earth Objects. As such, their survey cadence is built around optimizing NEO detection above all else. I'm interested in using their data to look for other objects also, but first I want to see how well HelioLinC can perform when tuned specifically for NEO detection. Can I recover the objects that CSS finds each night with HelioLinC? Does HelioLinC find anything more? Those questions are not answered in this study, but they're what's guiding the work long term. In this study I'll explore HelioLinC's capabilities with CSS' cadence using JPL's Small Body Identification API to pull data for NEOs only. As with my previous HelioLinC study, I'll add noise and fake detections to the data to see how HelioLinC responds to imperfect inputs.

I should emphasize that the results here won't apply to other surveys or telescopes. Survey cadence matters a lot and I'm tuning HelioLinC for NEO detection at CSS' very particular cadence. Anything I say about NEO detection really only applies to Catalina-like hardware and survey design with 4 observations over approximately 30 minutes. I would not expect the results of this study to apply to substantially different survey cadences or different telescopes.

Getting NEO observations with sbident

I wanted to pull in a universe of NEO observations and only NEO observations at CSS' cadence for testing. Once again, my sbident library is proving very useful. I used to just think of it as a tool for validating detections when I was using it with TESS data, but now I also think of it as a data source for generating synthetic observations for linking algorithms. Anyway, here's some example code that pulls NEO observations from JPL's Small Body Identification API at CSS' approximate cadence. I pulled a very wide field of view (25 degrees) just to get more NEOs into my dataset.

import numpy as np
from astropy.coordinates import SkyCoord
from sbident import SBIdent

tstart = 2458349.0 #starting jdate
times = np.linspace(tstart, tstart + 30.0/1440, 4) #4 observations over 30 minutes
center = SkyCoord(315, -15, unit="deg") #RA/DEC of center
mpc_code = '500' #geocentric observer MPC observatory code
hwidth = 25/2.0 #half width of field of view in degrees

for t in times:
  sbs = SBIdent(mpc_code, t, center, hwidth=hwidth, filters={'sb-group': 'neo'})
  #do stuff w/SSO observations contained in 'sbs.results'

As with my initial HelioLinC study, I will test both perfect observational data and observational data degraded by astrometric error and fake detections. The figure below shows perfect data on the left and imperfect data on the right.

Figure 1. The images above show the perfect (left) and imperfect (right) test data I use in this study. The left image shows JPL's RA and DEC positions for objects in a 25 degree field of view across 4 equally spaced observation times spanning 30 minutes. Because of the wide field of view and the short observation duration it's difficult to see the spatial separation in the observations, but they are separated. The image on the right is the same data except the observations have been corrupted by 1" σ noise and fake detections inserted at a multiple of 20x the number of true detections.

The NEO problem from a HelioLinC perspective

I'll start by describing what I think of as the dynamical characteristics of NEOs that differentiate them from other SSOs and which are pertinent to their detection with HelioLinC. This is my first foray into the study of NEOs; I'm still building my intuition for this class of solar system object.

First off, NEOs are defined as solar system objects with perihelion distances of less than 1.3 AU. While not always close to Earth, at some point in their orbit they do pass close to Earth's orbit. Due to their general proximity to Earth at the time of detection, NEOs tend to move relatively quickly on-sky for an Earth-based observer. This is useful because a NEO will traverse more pixels on a CCD for a given observation window than most other SSOs. As a result, a survey can observe a region of sky fairly briefly thus enabling more of the sky to be imaged per night. In a heliocentric frame this motion is not as pronounced, but NEOs are still closer than most SSOs to the Sun and as a result trace out relatively more of their orbit in time. These two features pair well with HelioLinC; you get comparatively more orbit arc out of a short duration observation window with NEOs than you would with other SSOs. Longer arcs generally yield better orbit determination and HelioLinC incorporates initial orbit determination in order to link detections.

Second, though related to the first point, because NEOs generally have larger eccentricities than main belt objects (Granvik et al. 2017, Malhotra, Wang 2016), their distance from the Sun (range) can change more quickly. In my previous HelioLinC study I was guessing constant, non time varying distances for the heliocentric range. With NEOs I suspect I may not be able to get away with that. By using a fixed range guess for all detections you're essentially circularizing your observations. This works surprisingly well a lot of times, but providing a range-rate guess (or guesses) to HelioLinC for NEOs might be necessary even with the small range variation implicit in the short observation window duration.

Results

With this rudimentary intuition of the NEO detection problem in place, I went about exploring the HelioLinC parameter space. While there are more parameters to tune at different levels of the HelioLinC algorithm, for this study I will focus on tuning the physical parameters of range and range-rate. To reiterate, range is the heliocentric distance of the object and range-rate is the change in the heliocentric distance of the object per unit of time.

To begin with, I tested a range guess of 1.3 AU and a range-rate guess of zero. Why not try the simplest thing first? For perfect astrometric data it works great. It's able to link all NEOs in the test data set with no confusions. However, it degrades quite rapidly when astrometric error is incorporated. At 1" σ of observational noise I was only able to link 21% (91/478) of the NEOs in the test data.

Note: there are 478 detectable objects in this test data set.

Table 1. Results: range=1.3 AU; range-rate=0 AU/day

Next I began experimenting with different range guesses. After trying many range guesses both larger and smaller than 1.3 AU, I found that guesses close to 1 AU performed quite well with both with perfect observational data and observations with astrometric error. However, I still couldn't get full recovery of the objects with the most noise added.

Table 2. Results: range=1.05 AU; range-rate=0 AU/day

Finally, I started pairing range guesses with non-zero range-rate guesses. After a few variations I found +/- 0.02 AU/day to yield good results. Specifically, I was guessing a range of 1.05 AU paired with range-rates of +0.02 AU/day, -0.02 AU/day and 0 AU/day. I ran these three range, range-rate pairs and then combined the links HelioLinC generated for each pair. With this I was able to recover at least 99.9% of the NEOs in all data I tested — including data with fake detections inserted.

Table 3. Results: range=1.05 AU; range-rates=[-0.2, 0, 0.2] AU/day

• Definitions:
• noise: the standard deviation of the normally distributed noise I'm adding to true object locations in arcsec.
• fakes: the number of fake detections inserted into the data expressed as a multiple of the number of true SSO sources in the observations.
• link: the unique set of detections in a cluster of HelioLinC propagated tracklets.
• HelioLinC links: the set of all individual HelioLinC generated links we're left with after removing subsets fully contained by another link.
• fully matched links: the count of links that contain all detections of an SSO. These are links where we were able to connect the same SSO detection to itself for every observation it was present.
• confusions: the count of links that contain more than one SSO or contain an SSO and one or more generated fake observations. This is a measure of error, the higher the count the worse it is.

The state space at the reference epoch

The core of the HelioLinC process is clustering detections that propagate to the same position and velocity at the reference epoch. Detections with the same position and velocity at the same reference time are the same object. The visual representation of this is shown below for the 3 range, range-rate pairs that achieve the best results. The final set of HelioLinC candidate links is the unique set of candidate links from each of these state space sets with subsets removed.

While the images may look like a sparse cloud of dots they're really a sparse cloud of dots plotted on top of one another with alpha blending darkening dots plotted on top of one another. These are the two detection tracklets propagating to the same position and velocity. If you click on an image (especially at range-rate=0 AU/day) to see a zoom you'll be able to find some instances of lighter dots where the tracklets did not propagate to the same position and velocity in significant quantity.

Figure 2. State space of common epoch propagated tracklets for range=1.05 AU; range-rate=-0.02 AU/day. The position component is on the left, the velocity component is on the right.

Figure 3. State space of common epoch propagated tracklets for range=1.05 AU; range-rate=0 AU/day. The position component is on the left, the velocity component is on the right.

Figure 4. State space of common epoch propagated tracklets for range=1.05 AU; range-rate=0.02 AU/day. The position component is on the left, the velocity component is on the right.

It's interesting how the state spaces with non-zero range-rate guesses have more widely dispersed velocities. This makes sense though. When you guess a range-rate of zero you're requiring the range to be the same at two different observation times. Circular orbits have a constant heliocentric range and thus the same velocity magnitude at all points along their trajectory. You can see that the velocity component of Figure 3 is clumped together in a fairly tight in a ball. The mean velocity of all of those dots is 28.9 km/s. That value isn't arbitrary; the mean velocity of Earth in orbit around the Sun is ~29.8 km/s. Since we were guessing a constant heliocentric range of 1.05 AU, it makes sense that the mean velocity of the tracklets we link would be slightly less than Earth's mean velocity.

Process improvements

Minimum link length constraint revised to minimum unique frame constraint

In my first HelioLinC study I was selecting clusters with at least N detections as candidate links. In this study I revised that so that I only select clusters with at least N distinct detection observation times as candidate links. This way I avoid generating candidate links that meet the minimum length threshold with multiple detections at the same observation time. For example, if my threshold was N=4 and there was a cluster with 4 detections in it but 2 of them were at the same observation time, I would not select that cluster as a candidate link.

No long links

I'm splitting long links (candidate links with more detections than the number of observation times) generated by the KD-tree with K-means clustering. I now run a loop on clusters with too many detections, splitting them with K-means until the longest length sub-cluster has no more than the number of observation time detections. I then select all sub clusters with at least N distinct detection observation times as above for candidate links. This way I never generate long links.

Tracklet pair generation

I generate tracklets by looking for pairs of observations within a certain span of one another. I wasn't really constraining that span intelligently before and as a result was testing way more tracklets than I needed to. I'm now constraining the span using a defined on-sky velocity and a duration that reflects the cadence of the observations. This yields a nice speedup by avoiding testing detection pairs that are silly.

6D clustering velocity normalization

I was previously using the maximum tracklet-joining velocity calculated by the solution to Lambert's problem to normalize all of the velocities for clustering. Oddly though, I found that I would get slightly better results when I ran test sets with fake data in them. Upon investigation I saw that the maximum velocity I was calculating for datasets with fake data inserted was always slightly higher than datasets without fakes. I'm specifying a maximum velocity above which I do not propagate tracklets back to the common epoch and both the fake and no-fake data had maximum velocities just below that (as they must given I was constraining them). But since they were slightly different they gave slightly different results. Now I'm just using that defined maximum velocity constraint to normalize all tracklet velocities prior to clustering.

Questions, ideas and next steps

Why is HelioLinC so effective with physical parameters that do not match those of the objects it links

This is still the fundamental open question for me. I don't mean that I expect HelioLinC not to work at all, just that it surprises me how well it works when we're choosing obviously 'wrong' physical parameters. For instance, not all of these NEOs are at a heliocentric range of 1.05 AU at the time of observation and it's highly unlikely that the 3 range-rate guesses closely match the true range-rates of any of the objects either. Still, HelioLinC manages to recover more than 99% of the NEOs for all of my test scenarios. I'm still curious to understand why HelioLinC is as robust as it appears to be.

Exploring the dynamical properties of objects that do and do not cluster

One obvious way to go about understanding HelioLinC's effectiveness would be to investigate the relationship of the guessed ranges and range-rates to the actual ranges and range-rates of the objects we're trying to link. The reason I haven't done this obvious thing yet is that I don't have the data. I should be able to pull this from JPL using one of their many useful APIs though.

Use real observation data from Catalina Sky Survey

Conversely, one way to disprove HelioLinC's effectiveness would be to use real data and find that it doesn't work as well as I think it does! There's always the possibility that the test data that I've constructed for this study doesn't really compare to real world data. In fact, I know my 'fake' detection insertion multiplier is just pulled out of thin air. I don't have an intuition for what the unique combination of hardware and image processing that CSS uses will yield in that respect. I'd like to soon be able to ingest their detection list and try HelioLinC on that instead of synthetic observations. It's very likely that real world data will be more challenging to deal with than naively modeled real world data.

Linking across fields in CSS data

Catalina Sky Survey visits a number of 'fields' each night and tries to link detections within a field to generate candidate links. As of now, CSS is not attempting to link across more than one field. Linking across fields may be something HelioLinC could enable if the combinatorial problem is not too substantial. By building up tracklets in adjacent fields you could perhaps link NEOs that were (for example) only detected in 2 frames of one field but then also detected in 2 frames of an adjacent field.

Conclusion

These case studies always prove so useful. They clarify my thinking and serve as helpful reference points on the problem landscape. What I initially considered a drawback (the short observation window duration) at the start of this study now seems like an advantage after testing NEO specific data. I was also surprised (at first) by how much range-rate guessing helped HelioLinC link things, but now I see how critical that will be for NEOs in particular. I think a NEO-tuned HelioLinC has real potential to be useful. Hopefully my next NEO study can incorporate more real data so I can do a comparison of what CSS finds versus a NEO-tuned HelioLinC.

Published: 3/27/2022