- Email: [email protected]
Automatic detection of geomagnetic sudden commencement via time–frequency clusters
Available online at www.sciencedirect.com
Advances in Space Research 48 (2011) 1537–1544 www.elsevier.com/locate/asr
Automatic detection of geomagnetic sudden commencement via time–frequency clusters Ali G. Hafez a,b,⇑, Essam Ghamry a,c a
National Research Institute of Astronomy and Geophysics (NRIAG), El-Marsad Street, Helwan, 11722 Cairo, Egypt b LTLab Inc., Kyushu University, 812-8581 Fukuoka-shi, Higashi-ku, Hakozaki 6-10-1, Japan c Space Weather Monitoring Center (SWMC), Helwan University, Ain Helwan 11795, Egypt Received 6 October 2010; received in revised form 15 May 2011; accepted 17 May 2011 Available online 26 May 2011
Abstract One of global processes in ionosphere–thermosphere–magnetosphere system is the geomagnetic storms. It is of great importance to develop an algorithm that auto-detects sudden commencement because it could be an indicator of onset of the geomagnetic storm. Automatic detection of geomagnetic sudden commencement is based on time–frequency clusters generated by spectrogram. Proposed algorithm is tested on data set collected from stations belong to the international real-time magnetic observatory network (INTERMAGNET). Maximum standard deviation of algorithm detection times is observed to be one minute of the corresponding arrival times published by National Geophysical Data Center (NGDC). Ó 2011 COSPAR. Published by Elsevier Ltd. All rights reserved. Keywords: Sudden commencement; Spectrogram; Automatic detection; SC timing; Wavelet; INTERMAGNET
1. Introduction Physical interaction of the interplanetary medium with Earth’s magnetosphere generates a variety of global processes in the ionosphere–thermosphere–magnetosphere system. One of these processes is the geomagnetic storm. A sudden enhancement of the solar wind dynamic pressure due to the interplanetary shock causes an abrupt increase of the H-component of geomagnetic field especially at low and middle latitudes (Araki et al., 1994). This phenomenon is called geomagnetic sudden commencement (SC). Geomagnetic storms are caused by the compression of strong electric field to Earth’s magnetosphere. This strong electric field is caused by a combination of solar wind velocity and southward interplanetary magnetic field, IMF Bz. The initial change of Earth’s magnetic field in this case is referred to storm sudden commencement (SSC).
⇑ Corresponding author at: LTLab Inc., Kyushu University, 812-8581 Fukuoka-shi, Higashi-ku, Hakozaki 6-10-1, Japan. E-mail address: [email protected] (A.G. Hafez).
After SSC occurs, the horizontal (H) component of Earth’s magnetic field in low latitudes increases like a step-function. This positive phase is called the initial phase and is followed by the main phase which is characterized by a depression in the horizontal (H) component lasting for several to 10 h. This depression caused by the westward direction ring current around Earth which can be measured by the Dst index. It is known that the SC is not a necessary condition for a storm to occur, and hence the initial phase is not an essential feature. But the essential feature of a storm is the significant development of a ring current. After reaching a minimum value, the horizontal component recovers slowly towards the initial undisturbed value. This phase is called the recovery phase. Geomagnetic storms usually last 1–2 days, but some may last for many days. High precision automatic detection of sudden commencement is of great importance, as geomagnetic storms could affect seriously the power grid like the power outage in Canada in 1989, radio and television interference and blackouts, hazards to orbiting astronauts and spacecrafts are caused by disturbances in the ionosphere and magnetosphere during geomagnetic storms.
0273-1177/$36.00 Ó 2011 COSPAR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.asr.2011.05.025
1538
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Automatic detection of onset timing of natural events is common in literature. Ali et al. (2009, 2010) used spectrogram and discrete wavelet transform (DWT) respectively to detect earthquake onset arrival timing. Mendes et al. (2005), detected variations of the horizontal component of the geomagnetic field related to geomagnetic storms by means of wavelets. Shinohara et al. (2005) developed an automatic real time detection and warning system of geomagnetic sudden commencements using real-time data from ground-based geomagnetic observations. They developed a method to detect the SC/SI (Sudden Impulse) using amplitude and increasing rate of geomagnetic field based on their statistical analysis of the SCs. Zossi de Artigas et al. (2008) applied daubechies function of the discrete wavelet transform to assess periodicities and long-term trends of geomagnetic storms from 1957 to 2004. Time–frequency representations (TFRs) such as spectrogram are two-dimensional tools for processing time-varying signals. Such signals occur in many real-world applications; examples include speech, music, biological signals, biomedical signals such as electrocardiogram (ECG) waveforms, impulse responses of wireless channels, radar, sonar acoustic waves and seismic waves. They are specifically designed to process time-varying signals as they jointly display time and frequency information demonstrating which frequencies occur at a certain time, or, at which times a certain frequency occurs. They combine both time and frequency domains information by displaying signals over a joint time–frequency (TF) plane. Through this paper we will use the spectrogram to detect SC. Time–frequency analysis establishes a local spectrum for any time instant, Allen and Rabiner (1977). High resolution cannot be obtained in both time and frequency domains simultaneously. The window must be chosen for locating sharp peaks or low frequency features because of the inverse relation between window length and the corresponding frequency bandwidth. Motivated by this fact, window length is chosen carefully in designing the algorithm for automatic detection of geomagnetic sudden commencement using frequency clusters generated by TFRs. In the proposed algorithm, different frequency clusters along time axis are determined. From these frequency clusters, third cluster is chosen to monitor the arrival of SC in it as the arrival is very clear in that cluster and threshold setting is easier as shown in Fig. 1. The recorded H-component of the geomagnetic field at KAK station is shown in (a). (b), (c) and (d) represent the first, second and third time–frequency clusters of this record respectively. 2. Time–frequency clusters One of the shortcomings of Fourier transform is that it does not give any information on time at which a frequency component occurs. Fast Fourier transform, FFT, analysis of non-stationary signals cannot describe the local transient features due to averaging over duration of the signal.
Fig. 1. Although SC is clear in all time–frequency clusters, but third cluster has the advantage of ease in determining a clear detection threshold level. (a) H component of the geomagnetic field at KAK station observed on January, 21, 2005. (b) First time–frequency cluster of (a). (c) Second time–frequency cluster of (a). (d) Third time–frequency cluster of (a).
That gave the motivation to have time frequency representations (TFRs). Processing in the time–frequency (TF) domain was first introduced with the introduction of Wigner distribution quadratic TFR (Wigner, 1932). In 1946, the sound spectrograph that relates to spectrogram TFR was introduced to represent speech signals for visual interpretation (Koenig, 1946). Cohen’s class of TFRs was provided based on Wigner distribution and its importance in signal processing (Cohen, 1995). Cohen’s class TFR is also defined in Allen and Rabiner (1977) and Altes (1980). Among all Cohen’s class TFRs and Wigner distribution, spectrograms are most often used for TF analysis due to their simplicity in theory and implementation. Spectrogram is the squared magnitude of the short time Fourier transform (STFT). STFT divides signal into successive short time segments by a sliding window and determines local spectrum by means of a Fourier transform of the segment (Zaman et al., 2003). It is based on the assumption that signal is stationary over short duration of the segment STFTfxðnÞg X ðm; xÞ ¼
1 X
xðnÞwðn mÞeixn ;
ð1Þ
n¼1
where m is the discrete-time index, x is the frequency parameter, x(n) is the signal to be analyzed, w(n m) is the windowed function centered at n = m. Matlab signal processing toolbox version 7 is used to determine STFT. The output STFT matrix has R rows
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
1539
and C columns, rows represent number of frequency clusters. Columns represent time points at which frequency clusters occur. R and C can be determined as follows:
time series, y(n), is divided to K segments, each segment contains eight samples, i.e.,
R ¼ bN =2c þ 1; C ¼ bðL NoÞ=ðW NoÞc;
for k = 1, . . . , K where K ¼ bL=8c. Then the STFT of each segment is:
ð2Þ
yk ¼ ðyð1 þ 8ðk 1ÞÞ; . . . ; yð8 þ 8ðk 1ÞÞÞT
bX c denotes the largest integer less than or equal to X : Y k ðxi ; mj Þ ¼ where N is the length of discrete Fourier transform, L and W are the length of the input time series and the analysis window respectively and No is number of samples the sections of the input time series overlap. Fig. 2 shows an SC recorded at AAE and the corresponding STFT. SC is clear in the two dimensional TF plane.
3. Methodology of the SC detector Duration of window width affects the TF localization offered by spectrogram. Furthermore, a trade-off exists between time and frequency localization because no window exists that is perfectly concentrated in both time and frequency domains due to uncertainty principle (Hlawatsch and Boudreaux-Bartels, 1992). A shorter analysis window yields better time localization but poorer frequency localization and vice versa. Time resolution is a very important aspect to make accurate SC detection therefore a shorter time window is used. In the proposed algorithm the input
8 X
yðn þ 8ðk 1ÞÞwðn mj Þe
ð3Þ
pffiffiffiffi 1xi n
;
ð4Þ
n¼1
where xi ¼ 2pði1Þ , mj = No(j 1), i = 1, 2, . . . , R and 8 j = 1, 2, . . . , C. The power spectrogram is: S k ðxi ; mj Þ ¼ jY k ðxi ; mj Þj
2
ð5Þ
Sk(xi, mj) is simply the squared value of STFT which denotes power of xith frequency cluster at time mj. Spectrogram is used in applications where high sampling frequency is used such as acoustic signals and detection of Pwave of earthquakes in Ali et al. (2009). The innovation in this research is to use the spectrogram to detect SC in geomagnetic data sampled at 1 sample per minute. The challenge in this paper is the setting of the spectrogram parameters to analyze such low sampling rate. The choice of these parameters is made to optimize between two aspects which are time and frequency resolution. The target is to use the shortest window which maintains good frequency resolution. Dividing the whole record to segments each one contains eight samples narrows the choices of N, W and No. N should have value not less than five to have good frequency
AAE-2005-Sept. 9 00.00-24.00 4
x 10 3.615
nT
3.61
3.605 3.6
3.595 2
4
6
8
10
12
14
16
18
20 3
7
2
Frequency
6 1 5 0
4 3
-1
2
-2
1
dB 0
5
10
15
20
Time (hours) Fig. 2. Time–frequency plane of a typical SC recorded at AAE station, upper panel is the H component of geomagnetic field at AAE station observed on September, 09, 2005. Lower panel is color map of the spectrogram of H component of that day. At abrupt changes higher power regions could be recognized.
1540
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Fig. 5. Amplitude of SCs of the data set used. Fig. 3. Detection threshold could be easily assigned to the third time– frequency cluster. (a-1) H component of the geomagnetic field at AAE station observed on March, 19, 2001. (a-2) SC of this event in time domain is zoomed in. (b-1) Third time–frequency cluster of (a-1). (b-2) SC of the third time–frequency cluster is zoomed in.
localization. Then values of W and No must be less than N, which is set to five. Consequently W and No are chosen to be four and three, respectively. Then we calculate the spectrogram matrix of every segment using following parameters, N = 5, Hanning window of width, W = 4, No = 3. The dimensions of this matrix, given by (2), are three rows and three columns. To reduce ripples which may cause false detection, average of each row will be calculated to form a column vector. Therefore each element of this vector corresponds to average power of a frequency cluster of the corresponding eight samples. After calculating spectrograms and their averages of all segments, we will have a matrix H(xi, k) which corresponds to the average powers of R frequency clusters for whole time series. These three clusters are shown in Fig. 1 in which we can recognize that it is difficult to use neither first nor second time–frequency clusters, (b) and (c)
Fig. 6. Period in minutes of the SCs used through this paper.
respectively, cause of difficulty in assigning threshold. Third frequency cluster, (d), is used to check the SCs of the data set
Fig. 4. Map of the stations used to collect data to test the proposed algorithm.
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Table 2 Original date and time of SCs used to compare proposed algorithm detection times with times made by Takano et al. algorithm.
50
Date
Time of SC
Shigeru et al. detection times
Proposed algorithm detection times
February 13 1990 February 15 1990 March 12 1990 March 20 1990 March 30 1990 April 09 1990 April 12 1990 April 17 1990 April 23 1990 May 18 1990
17:15
12:32
17:14
22:48
15:54 21:52 22:24 22:57
23:06
15:05 21:52 23:20 23:40
03:22
22:41
09:48
07:21 09:20 16:16
22:37
08:41 14:32 14:48 16:16 20:40 03:22 05:04 05:28 06:16 17:44 07:21
07:40
May 21 1990
10:22
May 26 1990
20:37
June 12 1990
08:20
July 28 1990
03:31
August 01 1990 August 26 1990 November 26 1990
07:41
17:15 01:28 21:57 22:35 15:04 15:03 22:44 12:59 07:21 09:30 08:43 16:27 05:06 03:45 07:20 19:33 05:27 05:07 13:04 13:36 10:21 22:15 20:37 23:25 20:59 22:58 03:31 07:20 07:41 14:11 11:05 14:07 23:32 03:01
nT/min
40
30
20
10
2.5
0
nT/min
20
40
60
80
100
Number of storms used Fig. 7. Maximum time variation of SCs nT/min of the data set used in this study.
Table 1 List of the observatories used in our study. Station name
Hartland Victoria Kakioka Lanzhou Qsaybeh Tamanrasset Phuthuy Guam Addis Ababa
1541
Code
HAD VIC KAK LZH QSB TAM PHU GUA AAE
Geographic
Geomagnetic
Lat. N
Long. E
Lat. N
Long. E
51.00 48.52 36.23 36.09 33.87 22.79 21.03 13.58 9.03
4.48 123.42 140.19 103.85 35.64 5.53 105.95 144.87 38.77
47.44 53.68 29.35 30.66 27.94 24.66 14.26 5.29 5.32
74.37 296.86 212.14 176.51 107.45 81.76 178.01 215.65 111.76
L-value
2.19 2.85 1.32 1.35 1.28 1.21 1.06 1.00 1.00
H ðxi ; kÞ ¼
06:26 15:03 22:44 07:21 08:43 03:26 07:19 03:36
05:43 23:32
23:35 01:58 21:51 22:40 23:21 22:14 22:43 00:26 07:20 00:36 10:30 13:48 05:23 05:30 00:42 08:46 05:25 19:40 07:40 22:13 20:45 23:03 20:38 23:27 18:05 14:59 05:58 06:31 07:42 07:45 05:43 11:50 10:22 19:58
K 1 X S k ðxi ; mj Þ K mj ¼1
04:38 21:55 03:37 07:39
03:38 05:28 06:56 11:20 15:12 07:38 08:16 15:28
02:03
22:33
21:56
20:33 20:56 21:52 22:08 23:28 00:58 02:58 08:18 08:58 14:58 01:06 03:30 06:00 06:40 06:56 07:38
20:58 05:56 11:22 06:12 01:57
05:46 06:08 07:12 08:10 08:56 23:32
ð6Þ
Fig. 3 shows averaged third time–frequency cluster together with H component of magnetic field as recorded at AAE station. Detection threshold is assigned as shown in figure which enables accurate SC detection. The computational time per one record (1 day contains 1440 samples) is less than one second.
4. Data set
Fig. 8. Probability distribution function (PDF) of the error in detection times generated by the algorithm comparing to detection times made by NGDC. Negative values means detection time made by algorithm leads detection time made by NGDC.
In this study, 111 sudden commencements followed by geomagnetic storms are checked. These events are between January 1993 and October 2008 based on the National Geophysical Data Center (NGDC) of the National Oceanic and Atmospheric Administration (NOAA). Checked data are 1 minute resolution of H-component (northward component) in low and middle latitudes from the international real-time magnetic observatory network
1542
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Fig. 9. Weak SC causes failure in detection. (a-1) H component of the geomagnetic field at KAK station observed on February, 15, 1990. (a-2) SC of this event in time domain is zoomed in. (b-1) Third time–frequency cluster of (a-1). (b-2) SC of the third time–frequency cluster is zoomed in. (c-1) H component of magnetic field at KAK station observed on May, 05, 1990. (c-2) SC of this event in time domain is zoomed in. (d-1) Third time–frequency cluster of (c-1). (d-2) SC of the third time–frequency cluster is zoomed in.
(INTERMAGNET). Fig. 4 shows location of all stations used in this paper. 5. Simulation results Data set collected from INTERMAGNET is chosen to satisfy different values of the following three parameters: 1- Amplitude of SC. 2- Period of SC increase. 3- Maximum time variation of SC. Shinohara et al. (2005) chose these parameters to set the limits of their real time SC detection algorithm to achieve a balance between SC/SI detection and the frequency of erroneous detection. They gave the following definitions of these parameters respectively: SC amplitude is the difference between onset magnetic field disturbance and the point at which SC is peaked. This peak point is usually exists 2–10 min after arrival of SC. Shinohara et al. (2005) specified 7.5 nT as the lower limit of SC amplitude for detection of SC. Applying our
proposed algorithm on the data set collected, it is found that our proposed algorithm could detect SC with magnitude 3.2 nT which is less than threshold set by Shinohara et al. (2005) which illustrates the capability of the proposed algorithm to detect small amplitude SCs as shown in Fig. 5. Period of SC increase is the time required (in minutes) for the SC to peak after initiation. Fig. 6 shows period of the SCs of data set under inspection. Maximum time variation of SC is the maximum change in magnetic field per minute observed within the period of SC increase. The smaller this variation rate is the harder to auto detect the SC time correctly. Shinohara et al. (2005) specified 2.5 nT/minute as the lower limit of maximum time variation of SC to enable detection of SCs. Fig. 7 shows that there are many events with slower rise time rate than 2.5 nT/min, minimum value was 0.98 nT/min. Such events with slow rise time could be detected by the proposed algorithm which means that our algorithm is very sensitive to very slow SCs. To calculate arrival time of SC, proposed algorithm in this paper is applied on 111 magnetic storms recorded by stations listed in Table 1. These detection times are
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
compared with times published on National Geophysical Data Center (NGDC). Number of false detections was 3 SCs which is a negligible number. Error is defined as the time difference between detection times calculated by proposed algorithm and those published by NGDC. Fig. 8 illustrates error probability distribution function. Average and standard deviation of error are found to be zero and 2 min, respectively.
6. Comparison study Proposed algorithm is tested on another data set used by Takano et al. (1999). Table 2 shows that proposed algorithm detection times are almost the same as (time of SC) which shows the efficiency of the proposed algorithm as it correctly detected almost all the storms based on simple algorithm using spectrogram except two SCs which could not be detected. Fig. 9 shows the two SCs which could not be detected by proposed algorithm. Failure of detection of these two events is cause of weakness of these SCs which make the amplitude at SC less than the detection threshold.
7. Conclusion A modern technique for determining geomagnetic storm sudden commencement is introduced using time– frequency representations. Third frequency cluster is used cause of its good synchronization timing with SC timing of original time series and ease of threshold setting. Simulation results elaborated that the algorithm is very sensitive to both small SC amplitudes and slow rise time rate. Applying the algorithm on data collected from INTERMAGNET, time of SC arrival computed by our proposed algorithm compared to arrival times published by NGDC, is observed to be insignificant. Proposed algorithm is applied on another data set used by Takano et al. (1999). Comparing results obtained from our proposed algorithm and results of Takano et al. (1999) proves that our proposed algorithm gave exact detection times although it is based on simple method using spectrogram instead of bi-orthogonal wavelet, which is more complicated.
Acknowledgements The results presented in this paper rely on data collected from magnetic observatories. We would like to thank the international real-time magnetic observatory network (INTERMAGNET) for promoting high standards of magnetic observatory practice (www.intermagnet.org). We are also indebted to Professor Shigeru Takano to provide us the data set they used in their study to achieve the comparison study in our paper.
1543
Appendix A This appendix contains the pseudo code of the proposed algorithm as follows: Read the H-component record Divide H-component record to segments Calculate number of segments Set the detection threshold While segment number is less than total number of segments Calculate third time–frequency cluster of the current segment If the value of the average value of the third cluster is larger than the detection threshold Calculate the timing of the segment Display that there is an SC in this record else Display that there is no SC in this record End if Add one to segment End while Read reference value of SC which is published by NGDC Define the error which is the difference between calculated SC timing and reference timing published by NGDC. References Ali, G. Hafez, Tahir, A. Khan, Kohda, T. Clear P-wave arrival of weak events and automatic onset determination using wavelet filter banks Digital Signal Process 20 (3), 715–723, 2010. Ali, G. Hafez, Tahir, A. Khan, Kohda, T. Earthquake onset detection using spectro-ratio on multi-threshold time–frequency sub-band. Digital Signal Process. 19 (1), 118–126, 2009. Allen, J.B., Rabiner, L.R. A unified approach to STFT analysis and synthesis. Proc. IEEE 65, 1558–1564, 1977. Altes, R.A. Detection, estimation, and classification with spectrograms. J. Acoust. Soc. Am. 67, 1232–1246, 1980. Araki, T. A physical model of the geomagnetic sudden commencement, in: Solar Wind Sources of Magnetospheric Ultra-Low-Frequency Waves, in: Engebretson, M.J. et al. (Eds.), . Geophys. Monogr., vol. 81. AGU, Washington, DC, pp. 183–200, 1994. Cohen, L. Time–Frequency Analysis. Prentice Hall, Englewood Cliffs, NJ, 1995. Hlawatsch, F., Boudreaux-Bartels, G.F. Linear and quadratic time frequency signal representations. IEEE Signal Process. Mag. 9, 1992. Koenig, W. The sound spectrograph. J. Acoust. Soc. Am. 18, 19–49, 1946. Zossi de Artigas, M., Fernandez de Campra, P., Zotto, E.M. Geomagnetic disturbances analysis using discrete wavelets. Geofı´s. Int. 47 (3), 257– 263, 2008. Mendes Jr., Margarete Oliveira Domingues, Domingues, Margarete Oliveira, da Costa, Aracy Mendes, Clua de Gonzalez, Alıcia L. Wavelet analysis applied to magnetograms: singularity detections related to geomagnetic storms. J. Atmos. Solar-Terr. Phys. 67, 1827– 1836, 2005. Zaman, Moushumi, Papandreou-Suppappola, Antonia, Spanias, Andreas. Advanced concepts in time–frequency signal processing made simple. In: 33rd ASEE/IEEE Frontiers in Education Conference, 2003.
1544
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Takano, Shigeru, Minamoto, Teruya, Arimura, Hiroki, Niijima, Koichi, Iyemori, Toshihiko, Araki, Tohru Automatic detection of geomagnetic sudden commencement using lifting wavelet filters. Lecture Notes in Computer Science, vol. 1721. Springer, Berlin, Heidelberg, 1999. Shinohara, M., Kikuchi, T., Nozaki, K. Automatic realtime detection of sudden commencement of geomagnetic storms. J. NICT 52 (3/4), 197– 205, 2005.
Solar-Geophysical Data, National Geophysical Data Center (NGDC), NOAA. Available from: http://sgd.ngdc.noaa.gov. Wigner, E.P. On the quantum correction for thermo-dynamic equilibrium. Phys. Rev. 40, 749–759, 1932.
Advances in Space Research 48 (2011) 1537–1544 www.elsevier.com/locate/asr
Automatic detection of geomagnetic sudden commencement via time–frequency clusters Ali G. Hafez a,b,⇑, Essam Ghamry a,c a
National Research Institute of Astronomy and Geophysics (NRIAG), El-Marsad Street, Helwan, 11722 Cairo, Egypt b LTLab Inc., Kyushu University, 812-8581 Fukuoka-shi, Higashi-ku, Hakozaki 6-10-1, Japan c Space Weather Monitoring Center (SWMC), Helwan University, Ain Helwan 11795, Egypt Received 6 October 2010; received in revised form 15 May 2011; accepted 17 May 2011 Available online 26 May 2011
Abstract One of global processes in ionosphere–thermosphere–magnetosphere system is the geomagnetic storms. It is of great importance to develop an algorithm that auto-detects sudden commencement because it could be an indicator of onset of the geomagnetic storm. Automatic detection of geomagnetic sudden commencement is based on time–frequency clusters generated by spectrogram. Proposed algorithm is tested on data set collected from stations belong to the international real-time magnetic observatory network (INTERMAGNET). Maximum standard deviation of algorithm detection times is observed to be one minute of the corresponding arrival times published by National Geophysical Data Center (NGDC). Ó 2011 COSPAR. Published by Elsevier Ltd. All rights reserved. Keywords: Sudden commencement; Spectrogram; Automatic detection; SC timing; Wavelet; INTERMAGNET
1. Introduction Physical interaction of the interplanetary medium with Earth’s magnetosphere generates a variety of global processes in the ionosphere–thermosphere–magnetosphere system. One of these processes is the geomagnetic storm. A sudden enhancement of the solar wind dynamic pressure due to the interplanetary shock causes an abrupt increase of the H-component of geomagnetic field especially at low and middle latitudes (Araki et al., 1994). This phenomenon is called geomagnetic sudden commencement (SC). Geomagnetic storms are caused by the compression of strong electric field to Earth’s magnetosphere. This strong electric field is caused by a combination of solar wind velocity and southward interplanetary magnetic field, IMF Bz. The initial change of Earth’s magnetic field in this case is referred to storm sudden commencement (SSC).
⇑ Corresponding author at: LTLab Inc., Kyushu University, 812-8581 Fukuoka-shi, Higashi-ku, Hakozaki 6-10-1, Japan. E-mail address: [email protected] (A.G. Hafez).
After SSC occurs, the horizontal (H) component of Earth’s magnetic field in low latitudes increases like a step-function. This positive phase is called the initial phase and is followed by the main phase which is characterized by a depression in the horizontal (H) component lasting for several to 10 h. This depression caused by the westward direction ring current around Earth which can be measured by the Dst index. It is known that the SC is not a necessary condition for a storm to occur, and hence the initial phase is not an essential feature. But the essential feature of a storm is the significant development of a ring current. After reaching a minimum value, the horizontal component recovers slowly towards the initial undisturbed value. This phase is called the recovery phase. Geomagnetic storms usually last 1–2 days, but some may last for many days. High precision automatic detection of sudden commencement is of great importance, as geomagnetic storms could affect seriously the power grid like the power outage in Canada in 1989, radio and television interference and blackouts, hazards to orbiting astronauts and spacecrafts are caused by disturbances in the ionosphere and magnetosphere during geomagnetic storms.
0273-1177/$36.00 Ó 2011 COSPAR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.asr.2011.05.025
1538
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Automatic detection of onset timing of natural events is common in literature. Ali et al. (2009, 2010) used spectrogram and discrete wavelet transform (DWT) respectively to detect earthquake onset arrival timing. Mendes et al. (2005), detected variations of the horizontal component of the geomagnetic field related to geomagnetic storms by means of wavelets. Shinohara et al. (2005) developed an automatic real time detection and warning system of geomagnetic sudden commencements using real-time data from ground-based geomagnetic observations. They developed a method to detect the SC/SI (Sudden Impulse) using amplitude and increasing rate of geomagnetic field based on their statistical analysis of the SCs. Zossi de Artigas et al. (2008) applied daubechies function of the discrete wavelet transform to assess periodicities and long-term trends of geomagnetic storms from 1957 to 2004. Time–frequency representations (TFRs) such as spectrogram are two-dimensional tools for processing time-varying signals. Such signals occur in many real-world applications; examples include speech, music, biological signals, biomedical signals such as electrocardiogram (ECG) waveforms, impulse responses of wireless channels, radar, sonar acoustic waves and seismic waves. They are specifically designed to process time-varying signals as they jointly display time and frequency information demonstrating which frequencies occur at a certain time, or, at which times a certain frequency occurs. They combine both time and frequency domains information by displaying signals over a joint time–frequency (TF) plane. Through this paper we will use the spectrogram to detect SC. Time–frequency analysis establishes a local spectrum for any time instant, Allen and Rabiner (1977). High resolution cannot be obtained in both time and frequency domains simultaneously. The window must be chosen for locating sharp peaks or low frequency features because of the inverse relation between window length and the corresponding frequency bandwidth. Motivated by this fact, window length is chosen carefully in designing the algorithm for automatic detection of geomagnetic sudden commencement using frequency clusters generated by TFRs. In the proposed algorithm, different frequency clusters along time axis are determined. From these frequency clusters, third cluster is chosen to monitor the arrival of SC in it as the arrival is very clear in that cluster and threshold setting is easier as shown in Fig. 1. The recorded H-component of the geomagnetic field at KAK station is shown in (a). (b), (c) and (d) represent the first, second and third time–frequency clusters of this record respectively. 2. Time–frequency clusters One of the shortcomings of Fourier transform is that it does not give any information on time at which a frequency component occurs. Fast Fourier transform, FFT, analysis of non-stationary signals cannot describe the local transient features due to averaging over duration of the signal.
Fig. 1. Although SC is clear in all time–frequency clusters, but third cluster has the advantage of ease in determining a clear detection threshold level. (a) H component of the geomagnetic field at KAK station observed on January, 21, 2005. (b) First time–frequency cluster of (a). (c) Second time–frequency cluster of (a). (d) Third time–frequency cluster of (a).
That gave the motivation to have time frequency representations (TFRs). Processing in the time–frequency (TF) domain was first introduced with the introduction of Wigner distribution quadratic TFR (Wigner, 1932). In 1946, the sound spectrograph that relates to spectrogram TFR was introduced to represent speech signals for visual interpretation (Koenig, 1946). Cohen’s class of TFRs was provided based on Wigner distribution and its importance in signal processing (Cohen, 1995). Cohen’s class TFR is also defined in Allen and Rabiner (1977) and Altes (1980). Among all Cohen’s class TFRs and Wigner distribution, spectrograms are most often used for TF analysis due to their simplicity in theory and implementation. Spectrogram is the squared magnitude of the short time Fourier transform (STFT). STFT divides signal into successive short time segments by a sliding window and determines local spectrum by means of a Fourier transform of the segment (Zaman et al., 2003). It is based on the assumption that signal is stationary over short duration of the segment STFTfxðnÞg X ðm; xÞ ¼
1 X
xðnÞwðn mÞeixn ;
ð1Þ
n¼1
where m is the discrete-time index, x is the frequency parameter, x(n) is the signal to be analyzed, w(n m) is the windowed function centered at n = m. Matlab signal processing toolbox version 7 is used to determine STFT. The output STFT matrix has R rows
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
1539
and C columns, rows represent number of frequency clusters. Columns represent time points at which frequency clusters occur. R and C can be determined as follows:
time series, y(n), is divided to K segments, each segment contains eight samples, i.e.,
R ¼ bN =2c þ 1; C ¼ bðL NoÞ=ðW NoÞc;
for k = 1, . . . , K where K ¼ bL=8c. Then the STFT of each segment is:
ð2Þ
yk ¼ ðyð1 þ 8ðk 1ÞÞ; . . . ; yð8 þ 8ðk 1ÞÞÞT
bX c denotes the largest integer less than or equal to X : Y k ðxi ; mj Þ ¼ where N is the length of discrete Fourier transform, L and W are the length of the input time series and the analysis window respectively and No is number of samples the sections of the input time series overlap. Fig. 2 shows an SC recorded at AAE and the corresponding STFT. SC is clear in the two dimensional TF plane.
3. Methodology of the SC detector Duration of window width affects the TF localization offered by spectrogram. Furthermore, a trade-off exists between time and frequency localization because no window exists that is perfectly concentrated in both time and frequency domains due to uncertainty principle (Hlawatsch and Boudreaux-Bartels, 1992). A shorter analysis window yields better time localization but poorer frequency localization and vice versa. Time resolution is a very important aspect to make accurate SC detection therefore a shorter time window is used. In the proposed algorithm the input
8 X
yðn þ 8ðk 1ÞÞwðn mj Þe
ð3Þ
pffiffiffiffi 1xi n
;
ð4Þ
n¼1
where xi ¼ 2pði1Þ , mj = No(j 1), i = 1, 2, . . . , R and 8 j = 1, 2, . . . , C. The power spectrogram is: S k ðxi ; mj Þ ¼ jY k ðxi ; mj Þj
2
ð5Þ
Sk(xi, mj) is simply the squared value of STFT which denotes power of xith frequency cluster at time mj. Spectrogram is used in applications where high sampling frequency is used such as acoustic signals and detection of Pwave of earthquakes in Ali et al. (2009). The innovation in this research is to use the spectrogram to detect SC in geomagnetic data sampled at 1 sample per minute. The challenge in this paper is the setting of the spectrogram parameters to analyze such low sampling rate. The choice of these parameters is made to optimize between two aspects which are time and frequency resolution. The target is to use the shortest window which maintains good frequency resolution. Dividing the whole record to segments each one contains eight samples narrows the choices of N, W and No. N should have value not less than five to have good frequency
AAE-2005-Sept. 9 00.00-24.00 4
x 10 3.615
nT
3.61
3.605 3.6
3.595 2
4
6
8
10
12
14
16
18
20 3
7
2
Frequency
6 1 5 0
4 3
-1
2
-2
1
dB 0
5
10
15
20
Time (hours) Fig. 2. Time–frequency plane of a typical SC recorded at AAE station, upper panel is the H component of geomagnetic field at AAE station observed on September, 09, 2005. Lower panel is color map of the spectrogram of H component of that day. At abrupt changes higher power regions could be recognized.
1540
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Fig. 5. Amplitude of SCs of the data set used. Fig. 3. Detection threshold could be easily assigned to the third time– frequency cluster. (a-1) H component of the geomagnetic field at AAE station observed on March, 19, 2001. (a-2) SC of this event in time domain is zoomed in. (b-1) Third time–frequency cluster of (a-1). (b-2) SC of the third time–frequency cluster is zoomed in.
localization. Then values of W and No must be less than N, which is set to five. Consequently W and No are chosen to be four and three, respectively. Then we calculate the spectrogram matrix of every segment using following parameters, N = 5, Hanning window of width, W = 4, No = 3. The dimensions of this matrix, given by (2), are three rows and three columns. To reduce ripples which may cause false detection, average of each row will be calculated to form a column vector. Therefore each element of this vector corresponds to average power of a frequency cluster of the corresponding eight samples. After calculating spectrograms and their averages of all segments, we will have a matrix H(xi, k) which corresponds to the average powers of R frequency clusters for whole time series. These three clusters are shown in Fig. 1 in which we can recognize that it is difficult to use neither first nor second time–frequency clusters, (b) and (c)
Fig. 6. Period in minutes of the SCs used through this paper.
respectively, cause of difficulty in assigning threshold. Third frequency cluster, (d), is used to check the SCs of the data set
Fig. 4. Map of the stations used to collect data to test the proposed algorithm.
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Table 2 Original date and time of SCs used to compare proposed algorithm detection times with times made by Takano et al. algorithm.
50
Date
Time of SC
Shigeru et al. detection times
Proposed algorithm detection times
February 13 1990 February 15 1990 March 12 1990 March 20 1990 March 30 1990 April 09 1990 April 12 1990 April 17 1990 April 23 1990 May 18 1990
17:15
12:32
17:14
22:48
15:54 21:52 22:24 22:57
23:06
15:05 21:52 23:20 23:40
03:22
22:41
09:48
07:21 09:20 16:16
22:37
08:41 14:32 14:48 16:16 20:40 03:22 05:04 05:28 06:16 17:44 07:21
07:40
May 21 1990
10:22
May 26 1990
20:37
June 12 1990
08:20
July 28 1990
03:31
August 01 1990 August 26 1990 November 26 1990
07:41
17:15 01:28 21:57 22:35 15:04 15:03 22:44 12:59 07:21 09:30 08:43 16:27 05:06 03:45 07:20 19:33 05:27 05:07 13:04 13:36 10:21 22:15 20:37 23:25 20:59 22:58 03:31 07:20 07:41 14:11 11:05 14:07 23:32 03:01
nT/min
40
30
20
10
2.5
0
nT/min
20
40
60
80
100
Number of storms used Fig. 7. Maximum time variation of SCs nT/min of the data set used in this study.
Table 1 List of the observatories used in our study. Station name
Hartland Victoria Kakioka Lanzhou Qsaybeh Tamanrasset Phuthuy Guam Addis Ababa
1541
Code
HAD VIC KAK LZH QSB TAM PHU GUA AAE
Geographic
Geomagnetic
Lat. N
Long. E
Lat. N
Long. E
51.00 48.52 36.23 36.09 33.87 22.79 21.03 13.58 9.03
4.48 123.42 140.19 103.85 35.64 5.53 105.95 144.87 38.77
47.44 53.68 29.35 30.66 27.94 24.66 14.26 5.29 5.32
74.37 296.86 212.14 176.51 107.45 81.76 178.01 215.65 111.76
L-value
2.19 2.85 1.32 1.35 1.28 1.21 1.06 1.00 1.00
H ðxi ; kÞ ¼
06:26 15:03 22:44 07:21 08:43 03:26 07:19 03:36
05:43 23:32
23:35 01:58 21:51 22:40 23:21 22:14 22:43 00:26 07:20 00:36 10:30 13:48 05:23 05:30 00:42 08:46 05:25 19:40 07:40 22:13 20:45 23:03 20:38 23:27 18:05 14:59 05:58 06:31 07:42 07:45 05:43 11:50 10:22 19:58
K 1 X S k ðxi ; mj Þ K mj ¼1
04:38 21:55 03:37 07:39
03:38 05:28 06:56 11:20 15:12 07:38 08:16 15:28
02:03
22:33
21:56
20:33 20:56 21:52 22:08 23:28 00:58 02:58 08:18 08:58 14:58 01:06 03:30 06:00 06:40 06:56 07:38
20:58 05:56 11:22 06:12 01:57
05:46 06:08 07:12 08:10 08:56 23:32
ð6Þ
Fig. 3 shows averaged third time–frequency cluster together with H component of magnetic field as recorded at AAE station. Detection threshold is assigned as shown in figure which enables accurate SC detection. The computational time per one record (1 day contains 1440 samples) is less than one second.
4. Data set
Fig. 8. Probability distribution function (PDF) of the error in detection times generated by the algorithm comparing to detection times made by NGDC. Negative values means detection time made by algorithm leads detection time made by NGDC.
In this study, 111 sudden commencements followed by geomagnetic storms are checked. These events are between January 1993 and October 2008 based on the National Geophysical Data Center (NGDC) of the National Oceanic and Atmospheric Administration (NOAA). Checked data are 1 minute resolution of H-component (northward component) in low and middle latitudes from the international real-time magnetic observatory network
1542
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Fig. 9. Weak SC causes failure in detection. (a-1) H component of the geomagnetic field at KAK station observed on February, 15, 1990. (a-2) SC of this event in time domain is zoomed in. (b-1) Third time–frequency cluster of (a-1). (b-2) SC of the third time–frequency cluster is zoomed in. (c-1) H component of magnetic field at KAK station observed on May, 05, 1990. (c-2) SC of this event in time domain is zoomed in. (d-1) Third time–frequency cluster of (c-1). (d-2) SC of the third time–frequency cluster is zoomed in.
(INTERMAGNET). Fig. 4 shows location of all stations used in this paper. 5. Simulation results Data set collected from INTERMAGNET is chosen to satisfy different values of the following three parameters: 1- Amplitude of SC. 2- Period of SC increase. 3- Maximum time variation of SC. Shinohara et al. (2005) chose these parameters to set the limits of their real time SC detection algorithm to achieve a balance between SC/SI detection and the frequency of erroneous detection. They gave the following definitions of these parameters respectively: SC amplitude is the difference between onset magnetic field disturbance and the point at which SC is peaked. This peak point is usually exists 2–10 min after arrival of SC. Shinohara et al. (2005) specified 7.5 nT as the lower limit of SC amplitude for detection of SC. Applying our
proposed algorithm on the data set collected, it is found that our proposed algorithm could detect SC with magnitude 3.2 nT which is less than threshold set by Shinohara et al. (2005) which illustrates the capability of the proposed algorithm to detect small amplitude SCs as shown in Fig. 5. Period of SC increase is the time required (in minutes) for the SC to peak after initiation. Fig. 6 shows period of the SCs of data set under inspection. Maximum time variation of SC is the maximum change in magnetic field per minute observed within the period of SC increase. The smaller this variation rate is the harder to auto detect the SC time correctly. Shinohara et al. (2005) specified 2.5 nT/minute as the lower limit of maximum time variation of SC to enable detection of SCs. Fig. 7 shows that there are many events with slower rise time rate than 2.5 nT/min, minimum value was 0.98 nT/min. Such events with slow rise time could be detected by the proposed algorithm which means that our algorithm is very sensitive to very slow SCs. To calculate arrival time of SC, proposed algorithm in this paper is applied on 111 magnetic storms recorded by stations listed in Table 1. These detection times are
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
compared with times published on National Geophysical Data Center (NGDC). Number of false detections was 3 SCs which is a negligible number. Error is defined as the time difference between detection times calculated by proposed algorithm and those published by NGDC. Fig. 8 illustrates error probability distribution function. Average and standard deviation of error are found to be zero and 2 min, respectively.
6. Comparison study Proposed algorithm is tested on another data set used by Takano et al. (1999). Table 2 shows that proposed algorithm detection times are almost the same as (time of SC) which shows the efficiency of the proposed algorithm as it correctly detected almost all the storms based on simple algorithm using spectrogram except two SCs which could not be detected. Fig. 9 shows the two SCs which could not be detected by proposed algorithm. Failure of detection of these two events is cause of weakness of these SCs which make the amplitude at SC less than the detection threshold.
7. Conclusion A modern technique for determining geomagnetic storm sudden commencement is introduced using time– frequency representations. Third frequency cluster is used cause of its good synchronization timing with SC timing of original time series and ease of threshold setting. Simulation results elaborated that the algorithm is very sensitive to both small SC amplitudes and slow rise time rate. Applying the algorithm on data collected from INTERMAGNET, time of SC arrival computed by our proposed algorithm compared to arrival times published by NGDC, is observed to be insignificant. Proposed algorithm is applied on another data set used by Takano et al. (1999). Comparing results obtained from our proposed algorithm and results of Takano et al. (1999) proves that our proposed algorithm gave exact detection times although it is based on simple method using spectrogram instead of bi-orthogonal wavelet, which is more complicated.
Acknowledgements The results presented in this paper rely on data collected from magnetic observatories. We would like to thank the international real-time magnetic observatory network (INTERMAGNET) for promoting high standards of magnetic observatory practice (www.intermagnet.org). We are also indebted to Professor Shigeru Takano to provide us the data set they used in their study to achieve the comparison study in our paper.
1543
Appendix A This appendix contains the pseudo code of the proposed algorithm as follows: Read the H-component record Divide H-component record to segments Calculate number of segments Set the detection threshold While segment number is less than total number of segments Calculate third time–frequency cluster of the current segment If the value of the average value of the third cluster is larger than the detection threshold Calculate the timing of the segment Display that there is an SC in this record else Display that there is no SC in this record End if Add one to segment End while Read reference value of SC which is published by NGDC Define the error which is the difference between calculated SC timing and reference timing published by NGDC. References Ali, G. Hafez, Tahir, A. Khan, Kohda, T. Clear P-wave arrival of weak events and automatic onset determination using wavelet filter banks Digital Signal Process 20 (3), 715–723, 2010. Ali, G. Hafez, Tahir, A. Khan, Kohda, T. Earthquake onset detection using spectro-ratio on multi-threshold time–frequency sub-band. Digital Signal Process. 19 (1), 118–126, 2009. Allen, J.B., Rabiner, L.R. A unified approach to STFT analysis and synthesis. Proc. IEEE 65, 1558–1564, 1977. Altes, R.A. Detection, estimation, and classification with spectrograms. J. Acoust. Soc. Am. 67, 1232–1246, 1980. Araki, T. A physical model of the geomagnetic sudden commencement, in: Solar Wind Sources of Magnetospheric Ultra-Low-Frequency Waves, in: Engebretson, M.J. et al. (Eds.), . Geophys. Monogr., vol. 81. AGU, Washington, DC, pp. 183–200, 1994. Cohen, L. Time–Frequency Analysis. Prentice Hall, Englewood Cliffs, NJ, 1995. Hlawatsch, F., Boudreaux-Bartels, G.F. Linear and quadratic time frequency signal representations. IEEE Signal Process. Mag. 9, 1992. Koenig, W. The sound spectrograph. J. Acoust. Soc. Am. 18, 19–49, 1946. Zossi de Artigas, M., Fernandez de Campra, P., Zotto, E.M. Geomagnetic disturbances analysis using discrete wavelets. Geofı´s. Int. 47 (3), 257– 263, 2008. Mendes Jr., Margarete Oliveira Domingues, Domingues, Margarete Oliveira, da Costa, Aracy Mendes, Clua de Gonzalez, Alıcia L. Wavelet analysis applied to magnetograms: singularity detections related to geomagnetic storms. J. Atmos. Solar-Terr. Phys. 67, 1827– 1836, 2005. Zaman, Moushumi, Papandreou-Suppappola, Antonia, Spanias, Andreas. Advanced concepts in time–frequency signal processing made simple. In: 33rd ASEE/IEEE Frontiers in Education Conference, 2003.
1544
A.G. Hafez, E. Ghamry / Advances in Space Research 48 (2011) 1537–1544
Takano, Shigeru, Minamoto, Teruya, Arimura, Hiroki, Niijima, Koichi, Iyemori, Toshihiko, Araki, Tohru Automatic detection of geomagnetic sudden commencement using lifting wavelet filters. Lecture Notes in Computer Science, vol. 1721. Springer, Berlin, Heidelberg, 1999. Shinohara, M., Kikuchi, T., Nozaki, K. Automatic realtime detection of sudden commencement of geomagnetic storms. J. NICT 52 (3/4), 197– 205, 2005.
Solar-Geophysical Data, National Geophysical Data Center (NGDC), NOAA. Available from: http://sgd.ngdc.noaa.gov. Wigner, E.P. On the quantum correction for thermo-dynamic equilibrium. Phys. Rev. 40, 749–759, 1932.