Tracking patchy pulsating aurora through all-sky images

Pulsating aurora is frequently observed in the evening and morning sector auroral oval. While the precipitating electrons span a wide range of energies, there is increasing evidence that the shape of pulsating auroral patches is controlled by structures in near-equatorial cold plasma; these patches appear to move with convection, for example. Given the tremendous and rapidly increasing amount of auroral image data from which the velocity of these patches can be inferred, it is timely to develop and implement techniques for the automatic identification of pulsating auroral patch events in these data and for the automatic determination of the velocity of individual patches from that data. As a first step towards this, we have implemented an automatic technique for determining patch velocities from sequences of images from the Time History of Events and Macroscale Interactions during Substorms (THEMIS) all-sky imager (ASI) and applied it to many pulsating aurora events. Here we demonstrate the use of this technique and present the initial results, including a comparison between ewograms (east– west keograms) and time series of patch position as determined by the algorithm. We discuss the implications of this technique for remote sensing convection in the inner magnetosphere.


Introduction
Diffuse aurora is the product of the pitch-angle scattering of magnetically bounce-trapped electrons and protons through wave-particle interactions (Davidson, 1990).While the most common and widely known aurora is the auroral arc, the diffuse aurora is perhaps more ubiquitous, and its most common type is pulsating aurora.Pulsating aurora is characterized by quasi-periodic variations in intensity and precipitating electrons with energies on the order of a few keV to several tens of keV (Johnstone, 1978).This type of aurora is most commonly seen in the morning sector auroral oval and persists for 1.5 h on average (Jones et al., 2011) but has been observed to last upwards of 15 h (Jones et al., 2013).Pulsating aurora often, but not exclusively, has an irregular patchy structure, an example of which can be seen in Fig. 1.The lifetime and size of patches is known to vary substantially, ranging from a few seconds to tens of minutes for the former and 10-200 km across for the latter (Royrvik and Davis, 1977).
Information about the nature of certain magnetospheric processes is reflected in the structure and evolution of the aurora they produce.For this reason, the auroral oval is often thought of as a screen onto which these processes project their dynamics, offering the opportunity to remotely sense activity in the magnetosphere through auroral observations, including magnetospheric convection.
Magnetospheric convection is the dominant mechanism of bulk plasma transport and circulation within the magnetosphere; it is driven by the coupling between the solar wind and Earth's magnetosphere (Axford, 1969;Angelopoulos et al., 1994).This bulk motion, also referred to simply as convection, is at its most basic level the result of E × B particle drift (Northrop, 1963) in the magnetosphere and ionosphere.It is often assumed that magnetospheric convection maps to the ionosphere due to magnetic field lines being treated as equipotentials under many circumstances.Potential drops parallel to the magnetic field will violate this assumption and introduce uncertainty into the mapping.
There are multiple ways of observing convection, including ground-based coherent and incoherent scatter radars as well as spacecraft.Coherent scatter radars such as the Super Dual Auroral Radar Network (SuperDARN) (Greenwald et al., 1978) measure the backscatter off of irregularities in Published by Copernicus Publications on behalf of the European Geosciences Union.ionospheric electron density, the motion of which is directly related to convection in the ionosphere.Meanwhile, incoherent scatter radars such as the Resolute Bay Incoherent Scatter Radar (RISR) (Doupnik et al., 1972) observe plasma motion through the detection of Thomson radiation.In situ measurements of ionospheric convection have been obtained from low-altitude satellites like the Defense Meteorological Satellite Program (DMSP) (Greenspan et al., 1986) and the Fast Auroral Snapshot Explorer (FAST) (Carlson, 1992).Higheraltitude spacecraft, including Cluster (Gustafsson et al., 1997) and the Time History of Events and Macroscale Interactions during Substorms (THEMIS) satellites (Angelopoulos, 2008), have measured magnetospheric convection.However, there are complications with all of these techniques.Su-perDARN benefits from extensive coverage, but as with all global observing systems, this comes at the expense of spatial and temporal resolution.Incoherent scatter radars have the opposite issue; they offer precise measurements but with a limited field of view (FoV).Satellites provide measurements of convection with high spatio-temporal resolution along their trajectories, but ultimately they are limited by their inability to separate variations in space and time.Convection is a challenging phenomenon to observe and additional techniques would be valuable.
There is increasing evidence that pulsating auroral patches are controlled by structures in the near-equatorial cold plasma (Rae, 2014).Since the motion of cold plasma is almost entirely determined by E × B drifting, these patches appear to move with ionospheric convection and could be used to create two-dimensional maps of convection (Yang et al., 2015).To create such maps, it is necessary to track large numbers of pulsating auroral patches.Past studies have tracked aurora by hand (e.g.Yang et al., 2015;Shiokawa et al., 2010), but there have been few attempts to apply an automatic algorithm.Blixt et al. (2006) developed a method using optical flow analysis that extracts object motion from sequences of images by assuming that variations in brightness are solely related to object motion.For pulsating aurora, this is a significant limitation due to their prominent pulsations and constantly evolving shape (Shiokawa et al., 2010).This study details a new method of routinely and quantitatively tracking auroral forms, focusing particularly on pulsating auroral patches.The technique provides information about persistent structures, including their geographic location and lifespan, located in sequences of all-sky images.

Dataset
The images of pulsating aurora processed for this paper were captured by the THEMIS all-sky imager (ASI) array (Donovan et al., 2006;Mende et al., 2008), the ground-based component of the NASA mission (Angelopoulos, 2008) designed to answer key questions about the aurora and substorms.Currently, the network consists of 21 ASIs stationed across northern North America, each capturing panchromatic, or "white light," images of the aurora on a 256 × 256-pixel CCD with a 3 s cadence.These instruments have been operating for over 10 years and have amassed tens of millions of images; on the order of 10 % contain pulsating aurora, which creates the ideal dataset for this study.
The rough estimate that 10 % of THEMIS images contain pulsating aurora can be calculated from the results of Jones et al. (2011), which state that up to 60 % of the images from the Gillam ASI contain pulsating aurora.Reducing this estimate by half to account for morning hours constituting half of the total operating hours of the camera and by half again to reflect cloud cover that is present approximately half of the time brings the estimate to 15 %.In addition, more poleward stations see pulsating aurora less frequently, thus decreasing the average to the order of 10 %.
The patches featured in this paper were imaged by the THEMIS ASIs located at Fort Smith, Sanikiluaq, and Whitehorse, Canada.Raw THEMIS images were used as inputs to the algorithm with data numbers cropped to 200 above the minimum value in the sequence being tracked and a maximum of 9000; these thresholds were determined through trial and error.The THEMIS ASI images and necessary programs for replicating this study are provided via Grono (2017).

Tracking algorithm
The tracking algorithm applied in this study is based on the work of Crocker and Grier (1996).Sequences of images are pre-processed through the application of a wavelet filter which convolves the images with a Gaussian surface.This operation transforms structures within the bandpass into Gaussian-like "blobs" and removes noise with a size outside of the bandpass, thus allowing auroral forms within the FoV of an ASI to be more easily identified by the tracking algorithm.The Crocker and Grier (1996) algorithm locates the peak intensity of each blob in x-y pixel coordinates by identifying the regional maxima which are within the upper 30th percentile of brightness across the entire image.This location is then refined to more accurately reflect the geometric centre of the structure by calculating the offset between the brightest pixel and the brightness-weighted centroid of neighbouring pixels.
Blobs in successive images are connected to form tracks under the assumption that they are non-interacting Brownian particles.To reduce the number of computations, a maximum distance is specified that is the farthest distance a blob is allowed to travel between frames.When the algorithm is attempting to identify where a blob moved, each candidate blob in the next image has a probability associated with it; this probability is the likelihood that a non-interacting Brownian particle would diffuse from the original location to the new location in the next image.By selecting the location with the greatest probability, blobs in successive images are clas- sified as the same structure.If a blob is not found in the next image -if it became too dim, for example -the algorithm keeps it in memory for a number of frames before declaring it to have permanently disappeared.This is an important feature for tracking pulsating aurora since without it the algorithm would tend to reidentify a patch as a new structure each time it leaves a dim state.
A number of algorithm parameters can be tuned, including the feature size that the filter enhances, the number of frames for which a blob can disappear or pixels it can travel before being identified as a different structure, and the minimum blob brightness and size.For the purposes of this study, the wavelet filter was set to enhance structures 11 pixels in diameter, remove noise on the order of 3 pixels, and perform median smoothing in a neighbourhood 5 pixels wide.The maximum distance blobs were permitted to travel between images was 10 pixels, in effect setting an upper limit to the patch velocities the algorithm could detect.For purely north-south motion relative to the FoV of the ASI, the limit was approximately 1700 m s −1 , while the limit for east-west motion was approximately 3700 m s −1 .Blobs were allowed to disappear for no more than five consecutive images; at the imaging cadence of the ASIs, this corresponds to 18 s, which is longer than the typical pulsation period of pulsating aurora.Figure 1 illustrates the effect of the filter as well as how blobs are identified.The minimum peak brightness of blobs was set to four; this number reflects the brightness of the byte-scaled image, not the raw ASI data.The minimum blob size is defined by the area under the blob and referred to as "mass", the minimum of which was set to 20 000 and is a reflection of the overall brightness of a feature.This parameter provides a threshold that helps to remove dim features in the camera which were not removed by the initial wavelet filtering.A final parameter informs the algorithm of the approximate diameter of the features of interest; this size was set to 17 pixels and provides an estimate of the typical size of pulsating auroral features observed by the cameras.

Results
Included in Fig. 2 are ewograms (east-west keograms) that demonstrate a difference between well-tracked (Fig. 2a) and poorly (Fig. 2b) tracked pulsating auroral patches.Each column in these figures was created by taking an ASI image and isolating the horizontal band of pixels between a patch's lifetime maximum and minimum latitude and flattening it down into a single row of pixels.The flattening was accomplished by manually choosing either the maximum or median value of each column depending on the brightness of the patch relative to its surroundings.When a patch was near another structure at least as bright as itself, the median value would often produce a superior ewogram by preventing the other structure from appearing in the ewogram instead of the desired patch.Chronologically arranging the columns for each image produces the final result.The white markers over-plotted onto the images indicate the longitude at which the tracking algorithm determined the patches to be located.It is important to note that the tracking algorithm can fail to identify a patch in an image, particularly when it becomes very dim, but reacquire it later; these instances are identified by gaps in the over-plotted markers.
The velocities shown in the top plots were calculated from the tracked auroral patches.To find the distance the patches travelled between any two frames, they were assumed to be at an altitude of 110 km and travelling at constant latitudes equal to the midpoint between the patches' true latitudes at subsequent time points.The time between frames was generally the 3 s resolution of the ASI, but it could be as high as 18 s if the patch was not identified by the algorithm for a maximum of five frames.Smoothing was necessary due to the presence of substantial jitter in the calculated patch locations.The smoothing technique utilized was a 30-data-pointwide moving boxcar average for which the uncertainties were estimated as the standard error of the speeds used to calculate the mean.
Moving beyond examples of isolated events, Figs.

Discussion
The results produced by this algorithm suggest that it can be well suited to tracking pulsating aurora.Given the correct parameters, it is capable of automatically tracking all patches within a sequence of images despite the pulsations characteristic of this type of aurora.In addition, patch speeds are low enough and the time resolution of THEMIS is high enough that the algorithm is generally able to accurately follow blobs between images.
While the algorithm can successfully track patches for extended periods if they maintain their shape or evolve slowly with time, structures that change dramatically or rapidly are problematic.Such patches can be recognized from their ewograms, the creation of which currently requires the patch to be tracked by hand with accelerations detected by the algorithm but not seen in the images.Figure 2b is an example of a patch with a shape that evolved substantially over the period it was tracked; notice that the white markers plotted on the ewogram suggest patch motion that is not seen in the ewogram itself.From the full ASI images for the lifetime of this patch, it was seen that the shape of the patch often evolved rapidly, affecting the location of the blobs produced by the wavelet filter.Thus, while the mismatch between the ewogram and the white markers in Fig. 2b suggests that the tracking algorithm did a poor job, the reality is that the variations in speed arose from the evolving shape.These types of patches are poorly suited for observing convection.
Conversely, the patch shown in Fig. 2a demonstrates that the algorithm is capable of accurately extracting the velocity of a patch from ASI images.Relative to the patch in Fig. 2b, that in Fig. 2a evolves quite slowly.A challenge to be overcome before pulsating aurora can be reliably used to observe convection is to separate slowly from rapidly evolving patches.
Moreover, it is possible for structure shapes to change so rapidly that they are quickly identified as a new blob; in fact, this is the most frequent outcome.As seen in Fig. 5, the most common length of time for a patch to be tracked was less than 20 frames.However, if the tracks of these short-lived patches are smooth enough, the data for these events have value in statistical studies as velocities can still be extracted from a nominally single patch which the algorithm may separate over time into two patches.In such a case, more than a single velocity would be calculated for a patch, but since it continually moves with the background convection speed, these velocities could still be used to observe ionospheric convection.
Regardless of how constant the shape of a patch remains, its pulsations introduce jitter into the output of the algorithm.Even minor variations in shape and brightness produce inconsistent results from the wavelet filter, creating fluctuations in the location of the blob peaks.In well-tracked patches, smoothing allows for velocities to be extracted which are comparable to what would be calculated by hand.Complicating matters further is that the elevation of a patch within the FoV of an ASI influences the quality of the track produced.This is due to the proportional relationship between the spatial resolution of the ASIs and the elevation within the images they produce; pixels at lower elevations capture aurora farther away than what is observed at zenith.Consequently, patches at lower elevations can exhibit unusual behaviour in one of two ways.A patch can either experience an extreme version of the jitter issue or have the opposite problem and not be observed to move whatsoever; both result from the large distances between pixel centres.Another issue with the output of the tracking algorithm at low elevations is that the wavelet filter always creates blobs along the edge of the FoV.These blobs have a consistent shape and can become the longest-tracked structures output by the algorithm.This issue arises because the night sky is brighter than the black pixels outside the FoV, and the filter enhances the edge of the FoV like it does to any of the auroral structures.This effect can be seen in Fig. 1.Fortunately, issues related to elevation can be resolved by restricting the patches studied to elevations greater than a threshold of approximately 50 • .
It must be mentioned that when multiple patches pass closely to one another or overlap, it is possible for the algorithm to switch from tracking one patch to another.Fortunately, this issue does not affect statistical studies or situations in which the velocity field is of interest.While these instances can affect a case study of a particular event, adjustments to certain algorithm parameters may remedy the issue.Patches coincident in time and nearby in space will not move with identical velocities, but it is often reasonable to assume that their motion will be similar.On average, neigh-bouring plasma blobs will convect at approximately equal velocities, but smaller-scale variations in plasma populations as well as electric and magnetic fields produce velocity vari-ations.Figures 3 and 4 demonstrate that the algorithm calculates patch velocities consistent with these expectations.On larger timescales, the velocities exhibit similar behaviour even though they may differ substantially at any particular moment.

Future work
The accurate tracking of the patches featured in Figs.2a, 3, and 4 promotes optimism that this algorithm will be successful if its complications can be mitigated or altogether avoided.Methods of improving the overall quality of the patch tracks that may warrant further exploration include aggressively limiting the results to higher elevations in the FoV of the ASIs, ignoring patches that were identified in too few frames for the data to be reasonably smoothed, and adjusting the input parameters of the bandpass filter to produce blobs which may be less sensitive to subtle changes in patch shape and brightness.Possibilities that have not yet been investigated include alternate filtering techniques and excluding patches with shapes that change too quickly based on the properties of their velocity time series.Instruments with improved spatial and temporal resolution would improve the tracks produced by this technique, chiefly by increasing the usefulness of tracks at lower elevations in the ASI FoV and decreasing the distance a patch is allowed to move between successive images.However, certain key areas would not benefit from this, including the jitter introduced by evolving patch structures and dim patches that can be temporarily lost by the algorithm.
Early results of patch filtration suggest that pulsating aurora can be separated into at least two distinct categories: one characterized by well-defined structures that follow convection and another more dynamic type with a motion that is difficult to ascertain and may be unrelated to structuring in cold equatorial plasma.This work is an initial step toward the ultimate goal of the automatic identification and tracking of the first type of pulsating aurora such that accurate maps of convection can be produced.Increasing the number of techniques available to observe and create maps of convection is beneficial to the community because it is a critical aspect of geospace dynamics which is difficult to observe.The coverage of the THEMIS ASI array is extensive.Under ideal meteorological and geomagnetic conditions, such as those described in Jones et al. (2013), it would be possible to produce a map of convection covering the majority of Canada and Alaska.Pulsating auroral patches arise from wave-particle interactions with cold plasma located in the region of the magnetosphere where magnetic field lines transition from tail-like to dipolar.When the next-generation auroral imaging network Transition Region Explorer (TREx) is deployed over the next few years, the combined data will offer an unprecedented view of this part of the magnetosphere.

Figure 1 .
Figure 1.A comparison of the output of the wavelet filter (a) to the original ASI image (b) in which blob locations are identified with green circles.The axes are in x-y pixel coordinates and the large purple circle identifies the patch that is the subject of Fig. 2a.

Figure 2 .
Figure 2. Smoothed longitudinal patch velocities (top) calculated from the locations over-plotted onto the bottom figures as white markers.The grey bars represent an estimated error in the speed, which was calculated by finding the standard error of the 30 values that went into the smoothed data point.Ewograms (bottom) comparing the actual patch motion to the result of the tracking algorithm, which is represented as white markers.
3 and 4 depict examples of multiple distinct pulsating auroral patches being tracked simultaneously within the FoV of an ASI.The former includes ewograms of three patches as well as a time series comparing the longitudinal component of their velocities.The latter figure considers velocity time series of three separate periods that include six or seven patches each; instead of ewograms, it includes ASI images containing markers denoting the locations of the majority of the patches.In both of the top two rows, a single patch is not marked in the ASI image since the lifetimes of the patches do not entirely overlap.Error bars are not shown in either of these two figures to avoid cluttering the images.

Figure 3 .
Figure 3.Time series (top image) of longitudinal velocity for three distinct, simultaneous patches that occurred between approximately 12:30 and 12:50 UT on 15 January 2011 at Fort Smith.Ewograms of the three patches (bottom three images).

Figure 4 .
Figure 4. Three time series (right) of longitudinal velocity for distinct, simultaneously tracked patches.Locations of the tracked patches (left) within the FoV of the ASI at an instant in time.In the top two rows, one event is missing from the ASI image due to the patches not being entirely simultaneous.

Figure 5 .
Figure 5. Histogram of patch lifetimes measured by the number of frames they appeared in for the 13:00 UT hour of THEMIS Whitehorse ASI on 28 December 2005.The bins are 20 frames wide.