Automatic detection of ionospheric Alfvén resonances using signal and image processing techniques

Induction coils permit the measurement of small and very rapid changes of the magnetic field. A new set of induction coils in the UK (atL = 3.2) record magnetic field changes over an effective frequency range of 0.1– 40 Hz, encompassing phenomena such as the Schumann resonances, magnetospheric pulsations and ionospheric Alfvén resonances (IARs). The IARs typically manifest themselves as a series of spectral resonance structures (SRSs) within the 1–10 Hz frequency range, usually appearing as fine bands or fringes in spectrogram plots and occurring almost daily during local night-time, disappearing during the daylight hours. The behaviour of the occurrence in frequency ( f ) and the difference in frequency between fringes ( 1f ) varies throughout the year. In order to quantify the daily, seasonal and annual changes of the SRSs, we developed a new method based on signal and image processing techniques to identify the fringes and to quantify the values of f , 1f and other relevant parameters in the data set. The technique is relatively robust to noise though requires tuning of threshold parameters. We analyse 18 months of induction coil data to demonstrate the utility of the method.


Introduction
Spectral resonance structures (SRSs) are a type of extremely low frequency magnetic field phenomena detectable in the 0.5-10 Hz region of the magnetic field spectrum.Their origin is attributed to the occurrence of ionospheric Alfvén resonances (IARs) between the E and F2 region layers of the upper atmosphere.The theoretical framework for IARs was originally posited by Polyakov and Rapoport (1981), and their existence was confirmed some years later by experimental detection in magnetic induction search coil data (e.g.Belyaev et al., 1989Belyaev et al., , 1990)).Further theoretical models have been developed (e.g.Lysak, 1993;Bösinger et al., 2009) to explain the phenomenon, and experimental work has provided evidence for their occurrence at low (Bösinger et al., 2002), middle (Kulak et al., 1999;Hebden et al., 2005) and high latitudes (Belyaev et al., 1999;Yahnin et al., 2003;Molchanov et al., 2004) in the Northern Hemisphere.Other researchers have shown SRSs are present in the Southern Hemisphere; for example, Yaun (2011) investigated SRSs measured at Halley Bay in Antarctica.
IARs are thought to arise from the partial reflection of magnetohydrodynamic Alfvén waves in the ionosphere, excited along magnetic field lines by the leakage of electric fields associated with terrestrial lightning activity (Greifinger and Greifinger, 1968;Fukunishi et al., 1996;Füllekrug et al., 1998).The waves reflect off the "walls" of the E region and F2 region where the gradients of electron density reach a maximum.This partial reflection sets up a series of resonant frequencies along the field lines, giving rise to between 2 and 12 harmonics (the so-called SRS).The SRSs are visible as a series of fringes in dynamic spectrograms of the magnetic field rate of change.
There are several parameters which can be observed from the SRSs which allow us to probe conditions within the ionosphere, but the primary one is the frequency interval ( f ) between spectral lines.This parameter can be used to determine length interval of the cavity if the corresponding foF2 values are available (e.g.Yahnin et al., 2003).From this other parameters such as the Alfvén velocity can be estimated.
The extraction of useful parameters (such as f ) from induction coil time series is time-consuming and has required a high level of manual intervention.Yahnin et al. (2003) used visual identification to produce values of occurrence rate and frequency interval or scale over a 3-year period from an instrument in Sodankylä, Finland.A combination of a manual point-and-click method and a semi-automatic method using fast Fourier transforms was applied by Hebden et al. (2005) to obtain key parameters, again from the Sodankylä data set.A paper by Odzimek et al. (2006) describes a semiautomated approach to capturing f using Bessel function fitting with a data set recorded in Poland.More recently, Yaun (2011) used a smoothing and sine curve fitting method to identify the values of the frequency interval in the aforementioned Sodankylä data set.However, none of the methods are fully automated, and all rely on some form of manual intervention.
In this study, we examine the use of a combination of signal and image processing techniques to identify the occurrence and properties of SRS fringes in spectral data.In Sect. 2 we describe the data used showing examples of the spectral features of interest, and detail the methodology used for automatic detection.Section 3 shows the results of automated processing of 18 months of data available from the UK observatory.We discuss the results in Sect. 4.

Data and methodology
In June 2012, the British Geological Survey (BGS) Geomagnetism team installed two horizontal high-frequency (100 Hz) induction coil magnetometers at the Eskdalemuir Observatory (55.3 • N, 3.2 • W; L = 3.2), in the Scottish Borders of the United Kingdom.The Eskdalemuir Observatory is one of the longest running geophysical sites in the UK (beginning operation in 1908) and is located in a rural valley with a generally quiet magnetic environment.
The instrumentation consists of two refurbished CM11-E induction coil magnetometers, north-south and east-west orientated, connected to a 24 bit Guralp digitiser.The digitiser converts the output signal from the coils for wired transmission to a logger located in a vault approximately 150 m from the coils.The data are recorded at 100 Hz by an onsite computer where they are collated into hourly files.The data are automatically collected from Eskdalemuir once per hour and permanently stored on the BGS network.Daily processing produces a set of spectrograph images for display on the BGS Geomagnetism website 1 .Since September 2012 the system has run almost continuously, losing only around 20 days of data (e.g.due to computer issues or component failure).

Visualising the SRSs
The SRSs appear as a set of fringes in time-frequency plots (i.e.spectrograms).Figure 1 shows three examples of spectrograms computed from the coil measurements which show 1 http://www.geomag.bgs.ac.uk/research/inductioncoils.htmlclear SRSs.To obtain these plots, the raw digitiser data are bandpass filtered between 0.5 and 10 Hz using a five-pole Butterworth filter in the time domain.This is primarily to remove the influence of longer period variations in the data, as the instrument level drifts with the general background magnetic field over a day.As the field is sampled at 100 Hz, there are no issues with aliasing or loss of power at these frequencies.
The raw data for 24 h are split into 862 slices of 100 s (10 000 samples) with an overlap of 100 s on either side.A Hanning window is applied to data to taper the edges to zero before the fast Fourier transform (FFT) is computed.The FFT data are scaled to SI units using the scaling and calibration factors of the digitiser and coil frequency response.
A spectrogram is created by stacking the outputs of 862 1-D FFT traces into a 2-D matrix.
The SRSs show a large variety of dynamic behaviour.On a typical day, they usually start as a single low frequency around the onset of local night-time before bifurcating into 5-10 separate fringes, increasing in frequency until around midnight (Fig. 1a).The fringes also widen in frequency until midnight before fading around the beginning of daylight.Occasionally, the fringes decrease in frequency slightly around 03:00 UT before eventually fading.No fringes are detected during daylight hours.
On rare occasions the SRSs can be seen at frequencies higher than the first Schumann resonance at ≈ 8 Hz (Fig. 1b).Other behaviours such as non-linearly spaced f occasionally appear (Fig. 1c).The spectrograms also show Pc1 pulsations (around 00:00-02:00 UT in Fig. 1b) and traces of horizontal broadband noise from "local" lightning storms within 2000 km (around 14:00 UT in Fig. 1c).The first Schumann resonance shows the typical diurnal variation related to the progressive triggering of tropical thunderstorms across the equatorial continental regions (i.e.Africa, America, Asia) as the Earth rotates.Despite being in an electrically and magnetically quiet region of the UK, there are still effects from harmonics of man-made noise, visible as thin vertical lines on the integer-valued frequencies (e.g. 1, 2, 3 Hz in Fig. 1b).

Method
We wish to automatically detect the position of the SRS fringes.The method can be divided into two parts.The first signal processing part approximately follows the methods of Odzimek et al. (2006) and Yaun (2011) to identify the SRS peak frequency within each FFT trace while the second part involves the use of image processing to link the individual peaks together and discern the continuous SRSs.The method is fully automatic once some thresholds have been set from inspection.As with any method, there is a trade-off to be made between detection of signal and sensitivity to noise.
We choose 14-15 February 2014 to demonstrate the method as it is a "good" day with low Kp and no local lightning activity or noise.The processing method starts with detrending and filtering the raw digitiser data from 0.5 to 10 Hz using a Butterworth five-pole filter in the time domain as shown in Fig. 2a.A series of FFTs are computed using the Welch method with 100 s of filtered data to produce 862 1-D spectra plots per day (from which the spectrograms are produced).The individual spectra are converted to SI units using the instrument response and digitiser calibration values (Fig. 2b).The smoothed non-stationary IAR peaks in each time slice are identified using the residuals from a best-fit sixth-order spline fit to remove the background trend (Fig. 2b).The size of the smoothing window depends on the number of samples per Hz.For this data set, a window of 14 samples (0.3 Hz) wide is sufficient to smooth the data.A peak-finding algorithm is used to identify the highest points in the detrended spectra (Fig. 2c).The peak-finding algorithm finds the maxima and minima of the data series and picks the peaks based upon the input threshold.The algorithm is set to ignore peaks where the ratio of the size difference between the largest and smallest peaks is less than 0.25, a value set after experimentation.Larger values imply the algorithm is more selective in finding peaks.The positions of the peaks from each time slice are then placed into a matrix which is treated as an image of "spots".In our processor, we use an FFT with a length of 4096 points for 100 Hz sample rate which gives a frequency resolution of 1/40.96Hz.The bandwidth of interest is between 0.5 and 10.5 Hz, meaning that the required number of points from each FFT is 410.Thus, the output "spots" image from the signal processing part has a size of 862 rows by 410 columns.This is used in the next stage which applies standard image processing techniques to link the "spots" together.
The next stage attempts to identify continuous SRS peaks by using standard image processing techniques (e.g.Kong and Rosenfeld, 1996;Gonzales and Woods, 2002).3a shows the input "spots" matrix to the image processing step.We to link together the continuous peaks but reject noise or sporadic peaks as best as possible.The "spots" images are first "dilated" by an extrusion of 15 pixels in the vertical direction and then "eroded" by a line reduction of 3 pixels in the horizontal direction.The two threshold values for the size of the dilation and the severity of the erosion were chosen from experimentation and relate mainly to the frequency and time resolution, as well as the general noise levels.Pixels that are touching after these operations are classed as been "connected" and are joined using the 2-D "bridge" morphological operator which dilates the connection.These regions are further dilated again to give the approximate position of the SRSs (Fig. 3b).
Next, the image in Fig. 3b and the original spectrogram are passed into the regionprops function in Matlab to determine the area and bounding box for each distinct region.The regionprops function works on binary images, splitting the connected positive-valued regions into distinct labelled components before computing various spatial properties such as the pixel area, the coordinates of the rectangular bounding box containing the region, the length of the perimeter and so on.
For this application, if the area of the detected region is less than 80 pixels, it is rejected and subsequently ignored.If the area is equal to or larger than 80 pixels, the bounding box is used to mask the spectrogram image and compute the weighted mean value of the pixels in that region of the spectrogram.This identifies the positions of pixels with the maximum amplitude in each row of the spectrogram matrix, and hence the peaks of the SRSs.These positions are shown in Fig. 3c.Values outside of local night-time (18:00-06:00) are ignored.Figure 3d shows the locations of the SRSs plotted onto the original spectrogram for the 14/15 February 2014.The positions of the SRSs can now be mapped as continuous lines throughout the spectrogram, and the relevant properties of interest can be derived from the data in the matrix shown in Fig. 3c.
The method is easy to implement and relatively quick to process, taking about 1 minute per day on a conventional desktop machine.In the next section, we show the results of applying the methodology to 18 months of induction coil data.

Results
Using the method described in Sect.2.2, we analysed the induction coil data set which consisted of 525 days (approximately 18 months) of induction coil data from September 2012 to February 2014.The data set is not perfect, as there are about 20 days in which data are completely missing.In addition, despite the site being relatively free of man-made interference, approximately 8 weeks show clear signs of contamination in the form of strong 1 Hz noise, which overwhelms the SRS features.To overcome these effects, the data set was split into "good" days by visual inspection to find 152 days where SRSs are clearly visible.As can be seen in Fig. 3d, the first Schumann resonance is also detected, so for our analysis, any SRSs detected above 6.5 Hz were ignored.
Figure 4 (upper) shows histograms of the number of fringes detected, the position in frequency in which the fringes occur and the f for all 525 days in the data set, while the lower panels show the results from the 152 selected "good" days.The histograms of the number of fringes suggest a modal value of 9 fringes for all the data, or 10 fringes for the selected data.Values of 1-2 fringes are most likely due to noise in the data sets.The frequency location of the fringes shows a preference between 0.5 and 4 Hz.SRSs occur slightly less often at frequencies of 4-6.5 Hz.In the upper middle panel there are larger values at the integer Hz frequencies.These are related to periods when man-made noise is prevalent in the data set.The right-hand histograms show the value of f across the data set.The modal value is between 0.5 and 0.75 Hz in both the full and partial selections of the data.This agrees well with previous observations from Yahnin et al. (2003), for example.
Figures 5 and 6 show the results of dividing the f data into the weighted average value per hour in UT for each of the 18 months.In Fig. 5 all the data were used, giving an equal number of data per month.In Fig. 6, the months are shown individually as each contains a different number of "good" days.The error bars show the 1σ variation of the data.
The average value of f varies seasonally.From Fig. 5, f is largest in Northern Hemisphere winter and smallest around the equinoxes.The error bars for the months of June and July are too large to make a definitive statement about the average value.Figure 6 shows a similar seasonal variation but suggests that f is smallest at the summer solstice.Note the end points of the analysis (e.g. at the start or finish of local night-time) have the largest error bars -particularly during Northern Hemisphere summer.
Although an obvious visual daily variation is the increase of f and f over the period from evening to midnight often followed by a decrease before daytime, this trend does not appear to be consistent enough to be observed in the hourly UTC averages.We would expect to see a convex-upward set of lines in Figs. 5 and 6.However, this is not observed, suggesting that it does not occur often enough to influence the overall statistics.

Discussion
Previous studies on the topic of SRS parameter evaluation have mostly focussed on deriving statistical properties from short periods or limited data (e.g.Bösinger et al., 2002;Hebden et al., 2005;Yaun, 2011) though Yahnin et al. (2003) had a data set covering over 4.5 years, Odzimek (2004) studied half a solar cycle, and Belyaev et al. (1997)  full solar cycle from 1985 to 1997.The variation of f both over diurnal and seasonal timescales, which we show in Figs. 5 and 6, has been noted by many other studies.Yahnin et al. (2003) also showed very strong indication of solar cycle variation, which we cannot yet validate.It is interesting to note that, in the Yahnin et al. (2003) data set, the frequency scale fell to its lowest value around 0.5 Hz (their Fig. 2b) at the peak of the solar cycle.The data presented here have been collected during solar maximum (albeit a weak cycle) and have a similar average frequency scale.However, there may also be a latitudinal effect as the Eskdalemuir site is located further south at L = 3.2 compared to L = 5.2.
As with any signal or image processing technique, a tradeoff exists between the sensitivity to signal and rejection of noise.The method is adept at detecting man-made interference, particularly linear features in the time-frequency plots.Figure 7 is an example of continuous noise detected on site for several weeks in May 2013, the cause of which is unknown.This type of interference accounts for the large valued bars in the integer frequency locations of the upper middle histogram in Fig. 4. Excluding such data is the easiest method for reducing their effect in statistical analyses.
Detection of SRSs is also sensitive to the threshold value set in the peak finding algorithm.Lowering the threshold increases the number of "spots" for the image processing step but introduces further noise and spurious SRSs. Figure 8 shows the effect of changing the peak-finding threshold from a high threshold value of 0.5 (Fig. 8a), where the ratio of the maximum to minimum peak range is 0.5, through the default value of 0.25 (Fig. 8b) to a low threshold value of 0.1 (Fig. 8c) for data from 13/14 June 2013.Visual inspection of the spectrogram image in Fig. 8a shows that there are some weak SRSs from about 22:00 to 00:00 (at 1-3 Hz).These are progressively picked out as the threshold value falls; however at threshold values of 0.1, most of the SRS features detected are just from noise.In general, such threshold parameters will have to be set by the user for individual data sets as they relate to the different instrument and processing settings.Although SRSs occur occasionally at frequencies higher than 7 Hz, it is difficult to distinguish when they overlap with the first Schumann resonance between 7.5 and 8.5 Hz.The analysis deliberately ignores SRSs which occur at frequencies higher than the first Schumann resonance, as the detected features there are generally due to noise.Finally, the method, as currently implemented, is imperfect at capturing fine SRS (e.g.Bösinger et al., 2004) where f < 0.2 Hz -though such features are rarely present in our data set.
We note that the method can be modified to detect other features of interest in data presented in the form of spectrogram by using a different bandpass filter and modifying the peak-finding threshold.For example, the method could be used to monitor the Schumann resonances or Pc1-type magnetosphere pulsations.The method can also be combined with other data sources such as the foF2 peak frequency from ionosonde data to infer near-real-time characteristics of the ionosphere.

Conclusions
We present an automatic method for detecting spectral resonance structures (SRSs) within the 0.5-6.5 Hz frequency range of induction coil data.The method is based on a mixture of signal and image processing techniques and can be used to locate SRS features and to quantify values such as the frequency interval ( f ) parameter in the data set.The technique is relatively robust to noise and works with low signal-to-noise ratio, though it does require user intervention to set a number of threshold values initially to suit the data set under investigation.We process 18 months of data to analyse the number of fringes, the occurrence location of SRSs in frequency and the frequency scale between fringes.These vary constantly throughout the year in a manner consistent with the results from previous studies.The method can also be applied to other phenomena of interest in the induction coil data such as detecting Schumann resonance properties or magnetospheric pulsations.

Figure 2 .
Figure 2. Signal processing protocol for identifying SRS peaks for 14-15 February 2014.See text for details.(a) De-trending and filtering of 100 s of 100 Hz raw digitiser data.(b) FFT of 100 s of Hanning-windowed data (blue), 14-point running average smoothed FFT data (red) and sixth-order spline background (black); (c) spline background removed from smoothed FFT with main peaks identified (circles).

Figure 3 .
Figure 3. Image processing protocol for identifying continuous SRS features for 14-15 February 2014.See text for details.(a) Image of peak "spots" identified in the smoothed FFT data; (b) dilation, erosion and bridging of "spots" image; (c) computation of weighted mean position of SRSs; (d) superposition of SRS peak locations with spectrogram.Features outside 18:00 and 06:00 are ignored in (b), (c), and (d).

Figure 4 .
Figure 4. Histograms of the number of fringes (left), frequency location (middle) and frequency interval ( f ) (right) for all 525 days (upper) and 152 selected "good" days (lower) from September 2012 to February 2014.Fringes with f > 6.5 Hz are ignored.

Figure 5 .
Figure 5. Plots of the hourly mean value of f across each month for all 525 days from October 2012 to February 2014.Data for October 2012-September 2013 (blue) and October 2013-February 2014 (red).Error bars show 1σ -level variation.

Figure 6 .
Figure 6.Plots of the hourly mean value of f for 152 selected "good" days in each month from September 2012 to February 2014.Error bars show 1σ -level variation.

FrequencyFigure 7 .
Figure 7. Example spectrogram of data from 25 to 26 May 2013 with very strong contamination from man-made interference (vertical lines).