By Paul Mohan

Abstract: This article examines the utility of applying image analysis methods to time-frequency spectrograms of RF meteor scatter. The constrained scatter-receiver geometry and associated Doppler behavior lends itself to efficient image processing approaches for extracting individual meteor scatter events from the raw, noise contaminated spectrogram. Classical segmentation and ‘blob’ analysis can be applied to derive metrics associated with event rates, Doppler ‘center of mass’ and proxies for scatter return energy.  A baseline feature extraction approach is described and applied to sample data collected at the Pine Hill Observatory (PHO) site in SE Michigan. Preliminary results presented here are hoped to encourage collaboration to explore potential correlations between spectrogram derived metrics and other more established measures used by the meteor observation community.

 

1 Introduction

Time-frequency spectrograms are an established tool supporting radio meteor observation, showing the dynamic time-dependent aspects of the scatter signal’s spectral structure. Real-time implementations find widespread use in audio analysis and are equally useful when applied to the down converted RF baseband signal from a meteor scatter receiver. Representative discussion and application of meteor observation using spectrograms are provided by Martinez (1988), Sanderson (2000), Bourdillon et al. (2005), Mohan (2003), and Verbelen (2020) to name a few. Figure 1 is an example of a meteor scatter spectrogram collected at PHO using illumination originating from the 55.154MHz video carrier of Canadian TV CHBX (46°35’39.84” N, 84°21’0.00” W) located 465km from our receiver site. The spectrogram’s horizontal axis spans 30 minutes (starting 9 January 2021; 16h00m UTC) and the vertical frequency axis spans 120Hz. The illumination carrier is centered on ~960Hz and vertical displacement of the scatter signals represent ± Doppler induced frequency shifts. Spectrograms typically use Fast Fourier Transforms with block lengths from 8192 to 32768 samples and sample rates from 8K to 48K samples/sec. As in all frequency domain analyses, there is a tradeoff between frequency and time resolution, with an improvement in one degrading the other.  The amplitude components displayed in the spectrogram are generally scaled to logarithmic power (dB) and can be color mapped using one of several common look-up tables. Several analysis programs are available providing spectrogram tools; Spectrum Lab, Spectran and R_Meteor to name a few.  It is upon these spectrogram ‘images’ that we investigate the application of image processing algorithms to enable extraction of scatter response metrics. Beyond event counting, we can define and quantify other metrics such as received scatter energy and Doppler ‘center of mass’.

 

Figure 1 – Time-frequency spectrogram from a 30 minute interval of moderate meteor activity observed at PHO.

 

There are several ongoing efforts in the radio meteor community seeking to exploit spectrograms and automatically extract potentially new and useful information. BRAMS is employing a neural network approach using training data manually extracted via crowd sourcing using spectrograms from continuous data collections by their network.

 

2 Approach

The functional flow of an algorithm developed at PHO for extracting meteor scatter metrics from their associated spectrograms is illustrated in Figure 2. Each of the algorithm steps is a standalone function generally available in image processing libraries with their foundational theory outlined in numerous texts such as Rosenfeld et al. (1976). While some of our algorithms internal parameter settings are driven by the particular collection system and signal characteristics at PHO, it is anticipated the basic structure of the algorithm is adaptable to other collection sites. Implementation of the processor was done using ImageJ, an open source and highly capable frame work / toolset for rapid prototyping of image analysis algorithms. Figure 3 shows the ImageJ processing script that executes the algorithm. It processes a typical spectrogram in 1 to 2 seconds on an i5 Core PC w/ 8GB RAM. Summary descriptions of the high-level processing steps follow.

 

Figure 2 – Functional processing flow of PHO algorithm for automated feature extraction.

Figure 3 – Image J processing script for automated extraction of scatter events.

Sub-image Extract

The raw spectrogram available to the processing chain generally contains other reference information not needed by the processor (e.g., axis labels, headers). Further, for the illumination source used by PHO, the transmitted video carrier contains 60 Hz frame rate harmonics which can cause vertically offset replications of strong scatter returns. The extraction of a fixed size sub-image, vertically centered on the zero-Doppler offset frequency, serves to ensure that only the desired window of the image is passed on for processing

Noise Filtering

RF noise sources and reflections from aircrafts occasionally pollute the meteor scatter spectrogram. While these are visually suppressed by the human observer, they must be removed to prevent false alarms within the downstream Segmentation stage. The approach chosen here is based on grayscale morphology, with Dougherty (1992) providing an introduction to its theoretical basis.  In our algorithm, we apply a vertical structuring element in what is termed an ‘Opening’ operation, which is erosion followed by dilation. This approach was driven by the relatively consistent structure of the Doppler signature of the desired scatter events which are strongly vertically oriented in the frequency dimension.  In contrast, the Doppler signatures from aircrafts are primarily horizontally oriented and extend over much longer time (horizontal) scales. These same features also apply to other forms of radio interference that produce horizontal signatures in the spectrogram. The Opening operation retains vertically oriented connected pixel groups within the spectrogram, and eliminates the remainder. It therefore also efficiently removes the salt and pepper contributions of random background noise.

Segmentation

Following extraction of appropriate vertical segments by the filtering stage, segmentation consolidates connected regions into single objects based on criteria of 1) minimum # of connected points in region and 2) amplitude threshold for inclusion in a region. In whole, this module serves to consolidate the associated vertical segments that make up the structures of interest in the spectrogram. This is particularly necessary for time-extended over dense trail returns that occupy larger areas within the spectrogram.

Thresholding

This step is simply a conversion from a grayscale image delivered by the segmentation stage into a binary image required by the Particle Analysis module to follow.

Particle Analysis

With the individual scatter regions built and isolated by segmentation, Particle Analysis then catalogs each region with a numerical identity and a bounding region descriptor that is placed into a region of interest (ROI) table. This list of regions is the key data set that defines the overlay map used in the following step.

ROI Overlay & Metric Generation

The ROI table generated in the previous step generates a region map that is overlaid back onto the raw sub-image extracted above. Figure 4 shows the output from this processing stage using the spectrogram of Figure 1 as input. Registration of the extracted regions (yellow outlines) with the original scatter traces is generally quite good, with the exception of the occasional false alarm. Within each overlaid region we can access the raw spectrogram data within it to measure specific metrics such as mean value, area, centroid, center of mass, among others. Table 1 provides an example of metrics generated for each overlay region. Here, X and Y are the image coordinates of the region centroid, XM, YM are coordinates of the center of mass and Width, Height are the dimensions of the region bounding box.

Figure 4 – Extracted and labeled scatter events overlaid onto the spectrogram of Figure 1.

 

Table 1 – Sample of region metrics associated with scatter events labeled in Figure 4.

Region Area Mean StdDev X Y XM YM BX BY Width Height
1 1112 77.44 42.18 107.49 187.51 107.58 197.32 104 0 8 320
2 306 58.04 28.25 1446.59 43.39 1446.6 43.13 1443 0 7 82
3 180 45.98 15.66 1437.4 44.9 1437.39 43.11 1436 7 3 80
4 30 44.1 25.1 1088.5 80 1088.5 76.75 1088 65 1 30
5 2970 79.37 41.15 1452.53 167.18 1452.33 174.23 1438 67 33 184
6 50 50.6 19.66 1503.98 84.48 1503.95 82.92 1503 71 2 26
7 383 64.13 29.51 1095.51 139.08 1095.63 140.4 1090 85 10 109
8 36 43.25 10.94 1437.5 108 1437.5 106.97 1437 90 1 36
9 57 40.96 14.48 1532.5 126.5 1532.5 122.85 1532 98 1 57
10 619 80.26 41.11 944.97 200.93 945.03 204.58 943 99 4 212
11 152 53.49 27.37 1609.72 131.12 1609.85 126.75 1608 100 3 65
12 394 54.35 23.03 1582.41 186.79 1582.47 189.58 1581 105 3 155
13 105 43.56 16.96 243.25 153.81 243.32 153.76 242 108 2 79
14 52 40.19 9.93 520 121 520.06 120.77 519 108 2 26
15 214 58.85 25.17 935.05 172.38 935.03 176.88 934 108 2 119

 

Post Analysis

Based on the ROI data extracted for each scatter event, various derivative analyses can be performed using tools like MatLab, Excel or MathCAD. Determining the appropriate interpretation and use of the data is subject to some considered thought accounting for the underlying scatter phenomenology and particular influences of the data acquisition chain. Some candidate analysis exhibits are presented below, recognizing they are placeholders as more data is analyzed.

 

3 Strawman Analysis Exhibits

Access to data extracted for each scatter event presents several options for further analysis. Initial application of the extraction algorithm has motivated several data reduction examples described below.

Received Scatter Energy

A metric representative of the scattered RF energy received per event can be considered based on the following:

  1. The spectrogram provides a measure of signal power (dB) mapped to 8-bits for display (e.g., 256 gray levels). The mean pixel amplitude within a given region is a proxy for power spectral density in units joules/(sec×Hz).
  2. The area of a given scatter ROI (in pixels) represents the product of frequency × time (Hz×sec)
  3. The product of (a) × (b) yields units of joules scaled by some constant k.

Applying this approach to the events extracted from Figure 1 yields the plot of Figure 5 showing the proxy for returned scatter energy versus time (30-minute interval). The abrupt cutoff on the low end of the energy scale is believed due to the limiting SNR of the receiver chain as well as effects of threshold settings currently used within the processing algorithm. Figure 6 shows the distribution of energy associated with the data of Figure 5. The presence of a double peak invites further investigation.

Figure 5 – Plot of received scatter energy (arbitrary units) versus time over the 30 min. interval of Figure 1.

 

Figure 6 – Distribution of received scatter energy (arbitrary units) over the 30 minute interval of Figure 5.

 

Event Doppler

The vertical displacement of the individual scatter returns in the spectrogram are a direct measure of relative motion induced Doppler shift. The zero Doppler reference line is centered on ~960 Hz and is the direct path signal (not scattered) from the illumination source. For each event region, the vertical pixel displacement (frequency offset) of its spectral ‘center of mass’ can be measured. Figure 7 shows a plot of Doppler versus time for the event regions indicated in Figure 4. A linear regression line shows the trend of Doppler frequency over time. Note that in the collection geometry available to PHO, the meteor trail segments being detected via reflection are nearly perpendicular to the horizontal line of sight to both the transmitter and receiver, in a configuration analyzed by Richardson (1999). Therefore, the radial motion of the scatter trail with respect to transmitter and receiver is relatively low, inducing only small Doppler frequency shifts (tens of Hz).

 

Figure 7 – Scatter event Doppler versus time.

 

Mean Power vs. Time-Bandwidth Product

A scatter plot of mean event received power versus the time-bandwidth product (pixel area) is illustrated in Figure 8. It reveals a trend (consistent in other spectrograms processed by PHO) that the larger the time-bandwidth product of an event, the greater the mean power. This would indicate that longer persisting events generally present a larger radar cross sections as well, versus a lower cross section that persists longer. The prominent clustering ‘wall’ observed at (time × bandwidth) values 25–28 is believed to be a lower bound imposed by the length of the structuring element used by the morphology filter (currently radius = 12 pixels), as well as the limiting frequency and time resolution currently in use. Reducing this effect is a future objective.

Figure 8 – Mean event power versus the time-bandwidth product (pixel area) over the 30 min. interval of Figure 1.

 

4 Comments and Next Directions

We have introduced examples of analysis templates and data extraction from event regions in meteor scatter spectrograms. A very limited set of spectrogram samples have been evaluated to date. Therefore, no general conclusions are promoted at this point. Moving forward, feedback from the broader meteor observation community would improve the value of these tools and methods. We anticipate that relationships could be drawn between the spectrogram derived metrics and those originating from other meteor analysis techniques. We further look forward to reporting results from spectrograms collected at PHO during the recent Quadrantid meteor shower. Preliminary processing shows solid detection performance despite dense populations of scatter events (~500 detections/30min during peak) as Figure 9 illustrates. Data reduction and analysis over extended time intervals is necessary to identify trends in behavior; for example, the anticipated bi-polar Doppler profile as the shower radiant rotates across the field of view.

Figure 9 – Awaiting further analysis , a continum of detections from the peak of the 2021 Quadrantid meteor shower.

Algorithm Performance:  Preliminary measurements of two key metrics of the event detection process, probability of detection (Pd) and false alarm rate (Pfa), yield values of Pd = 96%; Pfa = 5.6%.  These are based on a limited set of data and clearly need to be better assessed across a larger volume of spectrograms over varying environmental conditions and meteor activity levels. This is planned to be accomplished in the coming months.

We note that it is expected that different sets of threshold parameters within the algorithm could optimize detection/false alarm rates for a given class of meteor activity (i.e., sporadic events, meteor showers). Further, it is clear that the particulars of a specific data collection site and signal chain, (i.e., effective signal gain, time/frequency resolution, display intensity mapping, to name a few) will influence selection of algorithm parameters.

It is evident that on occasion aircrafts within the field of view which passing close to nadir produce reflections with relatively high Doppler slopes. This spectrogram signature mimics the vertical structure of meteor returns and false alarms are observed. Improvements in the filtering stage are needed to minimize these effects.

References

Bourdillon A., Haldoupis C., Hanuise C., Roux Y., Menard J. (2005). “Long duration meteor echoes characterized by Doppler spectrum bifurcation”. Geophysical Research Letters, 32. (5), L05805 (4 p.).

Dougherty E., (1992). “An Introduction to Morphological Image Processing”. SPIE Optical Engineering Press.

Martinez P. (1998). “Using Doppler DSP to study HF propagation”. RadCom, May 1998.

Mohan P. (2003). “Leonids 2002 high-resolution dynamic Doppler spectra using CW bistatic radar”. 2003 Leonid MAC Workshop, Poster session.

Richardson J. (1999). “The meteor meniscus: Meteor distance versus meteor zenith angle”. Meteor Trails, Journal of the AMS, (No. 4).

Rosenfeld A, Kak A., (1976). “Digital Picture Processing”, Academic Press.

Rueden C. T., Schindelin J., Hiner M. C. (2017). “ImageJ2 ImageJ for the next generation of scientific image data”. BMC Bioinformatics 18:529, PMID 29187165, doi:10.1186/s12859-017-1934-z (on Google Scholar).

Sanderson T. (2000). “Echoes from the Leonids”. RadCom, March 2000.

Verbelen F. (2020). “Radio meteors August 2020”. eMetN, 6, 73–80.