Wavelet Model of Geomagnetic Field Variations and Its Application to Detect Short-Period Geomagnetic Anomalies

Abstract

Geomagnetic data analysis is an important basis for the investigation of the processes in the near-Earth space, Earth magnetosphere, and ionosphere. The negative impact of geomagnetic anomalies on modern technical objects and human health determine the applied significance of the investigation and requires the creation of effective methods for timely detection of the anomalies. Priory complicated structure of geomagnetic data makes their formalization and analysis difficult. This paper proposes a wavelet model for geomagnetic field variations. It describes characteristic changes and anomalies of different amplitude and duration. Numerical realization of the model provides the possibility to apply it in online analysis. We describe the process of model identification and show its efficiency in the detection of sudden, short-period geomagnetic anomalies occurring before and during magnetic storms. Raw second data of the Paratunka and Magadan observatories and post-processed minute data were used in the paper. The question of noise effect on the proposed model results was under consideration.

1. Introduction

Geophysical monitoring data analysis is an important basis for investigating the processes in the near-Earth space (NES), Earth magnetosphere, and ionosphere. The recorded parameters of natural environments are used in many fundamental and applied investigations of solar-terrestrial physics and space weather. Of special actuality are the methods that provide online data processing and analysis and timely detection of anomalous phenomena in the NES. Anomaly detection and identification are crucial when constructing methods for the forecasting of natural catastrophic phenomena, such as magnetic storms, earthquakes, tsunamis, and powerful proton increases (GLE-events) etc.

The investigation object in the paper is the Earth magnetic field variations obtained by direct measurements at magnetic observatories of the international network (INTERMAGNET [1]). The collected sets of geomagnetic data allow one to investigate space-time features, study the dynamics of the processes in the Earth magnetosphere, and look for nonstationary manifestations of space weather. Thus, they are valuable satellite data support [2,3,4,5,6]. Investigation of geomagnetic data is of interest in different spheres of human activities, for example, research of atmosphere and hydrosphere [7], the study of the Earth deep structure (magnetotelluric sounding) [8,9], space weather forecast [2,3,4,5,6,10,11], seismic activity monitoring [12], and so on. The negative impact of geomagnetic disturbances on technical objects and on human health [13,14,15,16,17,18,19] determines the important applied significance of the research. The latest investigations show that sudden geomagnetic field short-period variations (geomagnetic pulsations) are a significant source of such an impact [11,12,20]. They occur in the polar region and propagate to the region of mid- and low-latitudes as geomagnetic activity grows. Thus, the task of the current diagnostics of such events is quite topical and is has been solved for several decades by different scientists and research groups [2,3,4,5,6,12,20,21,22,23,24].

The geomagnetic data modeling and analysis problems are associated with their complicated irregular structure and the presence of noises of different origin (natural and anthropogenic) [22,24,25,26]. The geomagnetic data structure and the anomalous feature occurrence times contain essential information on nonstationary processes in NES and characterize magnetic disturbance intensity at the place of data recording. Traditional methods for time series modeling (smoothing, spectral methods etc.) allow us to investigate long-term variations and regular diurnal changes of the magnetic field [2]. However, they are not effective enough to analyze sudden short-period features occurring during increased solar activity and magnetic storms [5,22,24,25,26]. At the present time, new methods and approaches are being developed to analyze complex geophysical data. For example, the authors of the paper [24] suggest using chirplet and warblet transforms to analyze fine structures of geomagnetic field variations. The authors [24] showed the possibility of estimating the instantaneous signal frequency and obtaining geomagnetic pulse characteristics based on a generalized wavelet transform. In the papers [5,6,25], it was proposed to apply a fuzzy logic technology to investigate geomagnetic data complicated structures. Application of fuzzy logic elements was further developed in the paper [22], where the authors consider a multi-scale approach and study morphological features of several intensive geomagnetic storms based on it. At the present time, to improve the efficiency of complicated raw data analysis, hybrid approaches, including machine learning methods, are frequently used [8,26,27,28]. For example, the paper [8] suggests a model of fuzzy wavelet neural network (FWNN), applying fuzzy logic elements and wavelet transform [29,30] to forecast a signal code transmitted from a well to the surface with wastewater. In the paper [28], the authors propose using deep learning structures to forecast the geoefficiency of coronal mass ejections (CME) and the times of geoeffective CME arrivals. They use a time series of satellite optical observations from 1996 to 2018 as the data for network training. However, the neural network apparatus validity and accuracy is determined by sampling representativity, and the required adaptation, taking into account geophysical data nonstationarity, decreases the neural network efficiency.

In this paper we propose to use the approach developed by the authors [31,32,33]. It is based on wavelet transform, which allows us to investigate complicated nonstationary changes in geophysical data and is widely used in physics, in particularly, in geophysics [4,7,9,26,34,35,36,37,38]. Based on the wavelet transform, we proposed an automated method to detect geomagnetic pulsations [34,35], constructed an algorithm for automatic detection of magnetic storms with increased risk of occurrence of geomagnetically-induced currents [3], and developed a method to calculate the geomagnetic activity index WISA [36,38]. Based on the discrete wavelet transform, the paper [4] proposed a method to determine the sudden beginning of a magnetic storm using the comparison of three methods for automatic detection of geomagnetic storms (first derivative analysis, Akaike information criterion, and discrete wavelet transform), the authors [3] show high efficiency of the discrete wavelet transform. The approach, suggested by the authors [3], using the algorithm of detecting coronal mass ejection shock fronts in ACE data on the solar wind before the trigger of storm arrival to the Earth, provides detection of a storm with high probability. The paper [7], applying statistical approaches and wavelet transforms, studies accelerometer measurements, solar index (F_10,7), geomagnetic index data (Kp), and different types of atmospheric density time variations. In the paper [26], the authors investigate two algorithms for detecting micro pulsations in geomagnetic data. The authors [26] showed that a combination of wavelet transform with deep neural network allows one to improve the detection efficiency of micro pulsations compared to discrete wavelet transform and threshold functions. In this paper we show the possibility to apply wavelet transform to optimize the calculation of Dst-index [33], which is a measure of field change caused by ring current occurring in the magnetosphere during magnetic storms [39]. The authors also developed a new technique for data wavelet decomposition. It allows one to suppress noise and to detect sudden short-period variations of the geomagnetic field and estimate their parameters [31]. This paper continues the investigations [31,32,33].

Compared to the traditional approaches, the proposed wavelet model of geomagnetic field variations allows us to describe regular changes and anomalous short-period variations of different forms and duration. They occur during increased geomagnetic activity and determine the field disturbance degree. The paper describes the model identification technique and shows the results of its application for data recorded at the Observatory Paratunka of IKIR FEB RAS (IAGA code is PET, 52.97N, 158.24 E, [1]). The method results are compared with the traditional approach based on Sq-variation application [40]. Sq-curve is used to describe the characteristic time variations of the geomagnetic field during quiet periods [40] and is an average smoothed curve obtained from the «quietest» field variations over a certain period (following the paper [40], five «quietest» variations over a current month are usually considered). Special attention in the paper is given to the possibility of detecting sudden short-period anomalous changes occurring at the preparation stage of a magnetic storm and may be considered as predictors. To estimate the method efficiency and study its capabilities in an online mode, the processing was performed for raw 2 Hz and one-minute geomagnetic data and for 2 Hz and one-minute data obtained after preliminary processing by magnetologists.

2. Description of the Method

In a wavelet space, geomagnetic field variations can be represented as

$$f\left(t\right)=Q^{char}\left(t\right)+D\left(t\right)+e\left(t\right)=Q^{char}\left(t\right)+{\displaystyle \sum }_{\rho }{\tau }_{\rho }^{pert}\left(t\right)+e\left(t\right),$$

(1)

where $Q^{char}\left(t\right)$ is a characteristic component describing geomagnetic field regular changes at an observation site (they are the changes of terrestrial magnetism elements with the period equal to the solar day duration [41]); $D\left(t\right)={{\displaystyle \sum }}_{\rho }{\tau }_{\rho }^{pert}\left(t\right)$ is the disturbed component including different-scale components of different form and duration and describing nonstationary anomalous changes occurring during increased geomagnetic activity (magnetic storm and substorm periods), $\rho$ is the component number; $e\left(t\right)$ is the noise component.

2.1. Identification of the Model Characteristic Component

In order to identify the model characteristic component $Q^{char}\left(t\right)$, we shall apply the multiple-scale analysis (MSA) [29,30,42]. We consider $L^2\left(R\right)$, the function space with finite energy (Lebesgue space [29,30]), as the signal space. Assume that the initial signal $f$ belongs to the subspace $V_j\subset L^2\left(R\right)$ of resolution $2^{-j}$:

$$V_j=clo{s}_{L^2\left(R\right)}\left({\varphi }_{j,n}:n\in Z\right), {\varphi }_{j,n}=2^{\frac{j}{2}}\varphi \left(2^{j}t-n\right), j\in Z.$$

Then, based on the MSA, we can obtain signal representation [29,30]

$$\begin{array}{l}{f}_j\left(t\right)=g_{j-1}\left(t\right)+g_{j-2}\left(t\right)+\dots +g_{j-m}\left(t\right)+f_{j-m}\left(t\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle {\displaystyle \sum }_{k=j-1}^{j-m}}{\displaystyle {\displaystyle \sum }_n}{d}_{k,n}{\Psi }_{k,n}\left(t\right)+{\displaystyle {\displaystyle \sum }_n}{c}_{j-m,n}{\varphi }_{j-m,n}\left(t\right).\end{array}$$

(2)

Each component in (2) is determined by the coefficient set ${\overline{c}}^j={\left\{c_{j,n}\right\}}_{n\in Z}\in V_j$, ${\overline{d}}^j={\left\{d_{j,n}\right\}}_{n\in Z}\in W_j$: $c_{j,n}=\langle f,{\varphi }_{j,n}\rangle$, $d_{j,n}=\langle f,{\Psi }_{j,n}\rangle$, $W_j$ is the wavelet space, ${\Psi }_{j,n}$ is the basic wavelet, $j$ is the scale, $m$ is the decomposition level. The coefficients ${\overline{c}}^j$ correspond to the signal approximating (smoothed) component $f_{j-m}$. The coefficients ${\overline{d}}^j$ determine the signal detailing components $g_{j-l}$.

Figure 1 shows the results of the MSA application to the one-minute magnetic data of the Paratunka observatory. The upper panel (Figure 1a–d) shows the results for 16 September 2021 (diurnal K-index = 5) and for 17 September 2021 (diurnal K-index = 25). The lower panel (Figure 1e–h) presents the results for the period of 1–30 September 2021. Diurnal values of K-index are indicated in the upper part of the panels. Smoothed components (Figure 1b,f) approximates signal trend and the detailing components (Figure 1c,d) contain short-period variations of different structures and duration. Analysis of the results shows that during calm periods (K-index has low values), MSA smoothed components allow us to obtain an approximation to the field calm variation. During disturbed periods, a significant increase of detailing component coefficients characterizes the field disturbance degree.

The smoothed component $f_{j-m}$ describes geomagnetic signal variations during quiet periods. The obtained approximation error depends on the decomposition level $m;$ the problem arises in determining the decomposition level $m^{*}$, which minimizes the error. Sq-curve can be used as a reference function, reflecting quiet geomagnetic field diurnal variation at an observation site [41]. As described above, Sq-curve is a mean smoothed curve of the five «quietest» intervals of the geomagnetic field over a current month. In this case, the approximation error in wavelet space can be determined as

$$U_m=\sqrt{{\displaystyle \sum }_{n=1}^L{\left|c_{j-m,n}-c_{j-m,n}^{Sq}\right|}^2},$$

(3)

where $c_{j-m,n}$ are the signal smoothed component coefficients at the decomposition level $m;$ $c_{j-m,n}^{Sq}$ are the Sq-curve smoothed component coefficients at the decomposition level $m;$ $n$ are readings; $L$ is the signal component length.

When choosing the decomposition level $m^{*}$, one should also take into account the fact that according to (2), during the change from the decomposition level $m$ to the decomposition level $m+1$, a part of the information on the signal turns to the detailing component $g_{j-l}$ (see (2)) that increases the losses

$$P_m=\sqrt{{\displaystyle \sum }_{n=1}^N{\left|c_{j,n}^m-c_{j,n}^{sig}\right|}^2},$$

(4)

where $c_{j,n}^m$ are the reconstructed signal coefficients obtained on the basis of wavelet recovery $f_{j-m}$ and using only approximation coefficients; $c_{j,n}^{sig}$ are the initial signal coefficients, $N$ is the signal length.

Thus, according to (3) and (4), we have a problem of multicriteria optimization in which the objective functions

$$\begin{gathered} {\Phi }^{U_m}={\displaystyle \sum }_{n=1}^N{\left|c_{j,n}^m-c_{j-m,n}^{m,Sq}\right|}^2\to min, \\ {\Phi }^{P_m}={\displaystyle \sum }_{n=1}^N{\left|c_{j,n}^m-c_{j,n}^{sig}\right|}^2\to min \end{gathered}$$

are mutually conflicting. In this case, the optimum level $m^{*}$ can be obtained by the method of constraint change ($\epsilon$-constraints):

$$\begin{gathered} {\Phi }^{U_m}={\displaystyle \sum }_{n=1}^N{\left|c_{j,n}^m-c_{j-m,n}^{m,Sq}\right|}^2\to min, \\ {\Phi }^{P_m^{\prime }}=\frac{1}{N}{\displaystyle \sum }_{n=1}^N{\left|c_{j,n}^m-c_{j,n}^{sig}\right|}^2\le \epsilon . \end{gathered}$$

(5)

The value $\epsilon$ is considered as an admissible level for ${\Phi }^{P_m^{\prime }}$, ${\epsilon }_i$ were estimated for each station $i$ as the dispersion of geomagnetic field variations during «quiet» periods (absence of magnetic storms and substorms).

From the above said we obtain the algorithm to choose the decomposition level $m^{*}$:

• We apply the MSA for the initial data estimated for each month of Sq-curves and obtain the representations $f_{j-m}\left(t\right)={{\displaystyle \sum }}_{n=1}^N{c}_{j-m,n}{\varphi }_{j-m,n}\left(t\right)$, $f_{j-m}^{Sq}\left(t\right)={{\displaystyle \sum }}_{n=1}^N{c}_{j-m,n}^{Sq}{\varphi }_{j-m,n}\left(t\right)$, $m=1,2,\dots ,J$, where $J\le {\log }_{2}N$ ($N$ is the signal length);
• For each decomposition level $m$ we carry out the reconstruction $f_{j-m}$ and $f_{j-m}^{Sq}:$$f_j^m\left(t\right)={{\displaystyle \sum }}_{n=1}^N{c}_{j,n}^m{\varphi }_{j,n}\left(t\right),f_j^{m,Sq}\left(t\right)={{\displaystyle \sum }}_{n=1}^N{c}_{j,n}^{m,Sq}{\varphi }_{j,n}\left(t\right)$ and estimate the error $U_m=\sqrt{{{\displaystyle \sum }}_{n=1}^N{\left|c_{j,n}^m-c_{j,n}^{m,Sq}\right|}^2}$ and losses $P_m=\sqrt{{{\displaystyle \sum }}_{n=1}^N{\left|c_{j,n}^m-c_{j,n}^{sig}\right|}^2}$;
• We determine the decomposition level $m^{*}$ providing the least error $U_m$ under admissible losses $P_m$ (conditions (5)).

The obtained estimates of $U_m$ and $P_m$ depend on the applied basic wavelet, the choice of which can be based on minimizing the loss function $P_m$ on the set of wavelet basis dictionaries. Figure 2 illustrates the calculation results of the losses for different basic wavelets. Daubechies wavelets of the 1-st–10-th orders, Db1–Db10 were used [30]. The estimates were calculated according to months (for January–June 2021) and for the whole period of the first half of 2021. To analyze geomagnetic activity, Figure 2 presents the K-index summary values for the periods under analysis. The results show close values of $P_m$ (see ratio (4)), including the decomposition 10-th level. A comparison with geomagnetic activity data shows loss increase when the geomagnetic field disturbance degree grows. The result indicates the complicated structure of geomagnetic data, containing features of different forms. Close loss values indicate the possibility of applying different wavelets for geomagnetic data approximation.

To assess the method efficiency, the results were compared with the method for construction of nonlinear approximations in wavelet basis described in the paper [42]. Decompositions into wavelet packages were applied. They allow one to obtain the best signal approximations of finite length by minimizing the concave cost function (entropy) [42].

When considering the dictionary $D$ of orthonormal bases

$$D={{\displaystyle \cup }}_{\lambda ϵ\mathsf{\Lambda }}{\mathfrak{B}}^{\lambda },$$

the approximation cost in the basis ${\mathfrak{B}}^{\lambda }={\left\{g_m^{\lambda }\right\}}_{1\le m\le N}$ can be estimated by Shur concave sum [43]

$$C\left(f,{\mathfrak{B}}^{\lambda }\right)={{\displaystyle \sum }}_{m=1}^N\Phi \left(\frac{{\left|\left(f,g_m^{\lambda }\right)\right|}^2}{\Vert f^2\Vert }\right), \Phi \left(x\right)=-x\ln x .$$

(6)

Then the optimal basis ${\mathfrak{B}}^{\alpha }$ minimizing the losses

$$P={\Vert f-\widehat{f}\Vert }_{L^2\left(R\right)’}$$

where $\widehat{f}$ is the estimate of the function $f$ is determined as

$$C\left(f,{\mathfrak{B}}^{\alpha }\right)=mi{n}_{\lambda ϵ\Lambda }C\left(f,{\mathfrak{B}}^{\lambda }\right).$$

In the package space [42]

$$W_j^p=W_{j+1}^{2p}\oplus W_{j+1}^{2p+1}$$

the optimal basis is [43]

$$W_j^p=\left\{\begin{array}{c}{W}_{j+1}^{2p}\cup W_{j+1}^{2p+1}, if C\left(f,W_{j+1}^{2p}\right)+C\left(f,W_{j+1}^{2p+1}\right)<C\left(f,{\mathfrak{B}}_j^p\right),\\ {\mathrm{W}}_j^p, if C\left(f,W_{j+1}^{2p}\right)+C\left(f,W_{j+1}^{2p+1}\right)\ge C\left(f,{\mathfrak{B}}_j^p\right).\end{array}\right.$$

(7)

Recursive computation of bases (7), when moving upwards along a wavelet package tree, allows us to obtain an optimal basis, minimizing the cost (6). Construction of an optimal basis is described more in detail in this paper [43].

Estimate results, obtained according to the data of the Paratunka observatory for January 2021 (K-index total value = 196) and for March 2021 (K-index total value = 371), are illustrated in

Table 1. Method efficiency estimate results, obtained from the data for January 2021.

Wavelet	$\mathbf{Losses} {\mathit{P}}_{\mathit{m}} \mathbf{of} \mathbf{the} \mathbf{Suggested} \mathbf{Method}$	$\mathbf{Losses} \mathit{P} \mathbf{of} \mathbf{the} \mathbf{Optimal} \mathbf{Basis}$	$\mathbf{Errors} {\mathit{U}}_{\mathit{m}} \mathbf{of} \mathbf{the} \mathbf{Suggested} \mathbf{Method}$	$\mathbf{Errors} \mathit{U}$ of the Optimal Basis
Db1	396.09	165.05	21.05	74.06
Db2	349.63	141.32	21.37	53.99
Db3	331.57	131.04	20.99	51.76
Db4	330.26	130.85	21.26	50.54
Db5	329.77	130.82	21.31	50.80
Db6	327.42	130.91	21.06	51.55
Db7	326.31	130.87	21.22	51.28
Db8	329.35	131.02	21.32	50.06
Db9	324.30	130.90	21.14	50.89
Db10	324.20	130.96	21.17	51.20

and

Table 2. Method efficiency estimate results, obtained from the data for March 2021.

Wavelet	$\mathbf{Losses} {\mathit{P}}_{\mathit{m}} \mathbf{of} \mathbf{the} \mathbf{Suggested} \mathbf{Method}$	$\mathbf{Losses} \mathit{P} \mathbf{of} \mathbf{the} \mathbf{Optimal} \mathbf{Basis}$	$\mathbf{Errors} {\mathit{U}}_{\mathit{m}} \mathbf{of} \mathbf{the} \mathbf{Suggested} \mathbf{Method}$	$\mathbf{Errors} \mathit{U}$ of the Optimal Basis
Db1	495.86	185.22	21.52	126.27
Db2	459.81	170.22	21.34	68.39
Db3	439.88	156.99	21.52	74.70
Db4	438.27	155.61	21.57	63.92
Db5	437.79	148.49	21.39	62.56
Db6	437.09	146.00	21.51	61.23
Db7	436.67	142.27	21.57	53.66
Db8	336.82	145.82	21.41	61.57
Db9	437.02	145.80	21.44	61.90
Db10	436.73	140.62	21.57	52.87

, respectively. They show that the losses $P$ of the optimal basis are significantly less than the losses $P_m$ of the suggested method. However, the errors $U$ of the optimal basis determined according to ratio (3) as

$$U={\Vert \widehat{f}-f^{sq}\Vert }_{L^2\left(R\right)’}$$

significantly exceed the errors $U_m$ of the suggested method, where $f^{sq}$ is the Sq-curve; as long as the Sq-curve, which describes the calm course of geomagnetic field variations, is taken as a reference function for the main condition in method (5) when making estimates and is the one providing the least error $U_m$. In this case, the decomposition level $m$ is determined so that the losses $P_m$ do not exceed the acceptable value $\epsilon$ (see (5)). The result confirms the efficiency of the suggested method.

Comparison of the results obtained for different basic wavelets (Db1–Db10) shows close values of errors $U_m$ and losses $P_m$ both during low (January 2021, Table 1) and during high (March 2021, Table 2) geomagnetic activity. The results confirm the author’s statement [42] that if a signal contains different structure features localized at different times, it is impossible to construct a basis adapted to all the structures. That also confirms the possibility of application of different wavelets for geomagnetic data approximation.

2.2. Identification of the Model Disturbed Component

Geomagnetic field disturbance degree is the quantity that, according to the paper [41], can be estimated by finding the difference between the highest and the least deviations of signal values from the corresponding values of a quiet diurnal variation. In the proposed model (1), geomagnetic disturbances are described by the component $D\left(t\right)={{\displaystyle \sum }}_{\rho }{\tau }_{\rho }{}^{pert}\left(t\right)$. Following the paper [31] and considering the significant nonstationarity $D\left(t\right)$, we shall use the continuous wavelet transform and adaptive threshold functions to identify it.

Continuous wavelet transform is determined as [29,30]

$$W_{b,s}={\left|s\right|}^{\frac{-1}{2}}{{\displaystyle \int }}_{-\infty }^{+\infty }f\left(t\right)\Psi \left(\frac{t-b}{s}\right)dt;s,b\in R,s\ne 0.$$

Then for $D\left(t\right)$ we obtain the representation:

$$D\left(t\right)={{\displaystyle \sum }}_{b,s}{E}_{1,s}\left(W_{b,s}\right){\Psi }_{b,s}\left(t\right)+{{\displaystyle \sum }}_{b,s}{E}_{2,s}\left(W_{b,s}\right){\Psi }_{b,s}\left(t\right),$$

$E_{1,s}\left(\cdot \right), E_{2,s}\left(\cdot \right)$ are the threshold functions determining the respectively weak (of $i=1$ class) and strong (of $i=2$ class) geomagnetic disturbances

$$E_{1,s}\left(W_{b,s}\right)=\left\{\begin{array}{c}0, if \left|W_{b,s}\right|\le T_{1,s} or \left|W_{b,s}\right|>T_{2,s}\\ W_{b,s}, if T_{1,s}<\left|W_{b,s}\right|\le T_{2,s}\end{array}\right.;$$

(8)

$$E_{2,s}\left(W_{b,s}\right)=\left\{\begin{array}{c}0, if \left|W_{b,s}\right|\le T_{2,s} \\ W_{b,s}, if \left|W_{b,s}\right|>T_{2,s}\end{array}\right.,$$

$T_{i,s}=U_i\times S{t}_s$, $i=1,2$, $U_i$ is the threshold coefficient, $S{t}_s=\sqrt{\frac{1}{l-1}{{\displaystyle \sum }}_k^l{\left(W_{b,s}-\overline{W_{b,s}}\right)}^2}$, $\overline{W_{b,s}}$ is the wavelet coefficient average calculated in a moving time window of length $l$.

Taking into account the equivalence of the discrete and continuous wavelet transform [29,30], when selecting the wavelet $\Psi$ we can use the estimates obtained for discrete decompositions (Table 1 and Table 2, Section 2.1).

The intensity of positive and negative disturbances of the geomagnetic field of each class $i$ at the time $t=k$ can be estimated as

$$I_k{}^{i\pm }\left(k\right)={\displaystyle \sum }_s\left|W_{b,s}^{\pm }\right|.$$

(9)

The thresholds $T_{i,s}$ in (8) are determined by the root-mean-square deviation ($S{t}_s$), which is estimated dynamically in a moving time window of duration $l$. That provides adaptation of thresholds $T_{i,s}$ to signal properties change.

The coefficients of the thresholds $U_i$ (see (8)) can be determined by minimizing the posterior risk [44]. In that case, having posterior data losses can be calculated as

$$R_s\left(y\right)={\displaystyle \sum }_{z=0}^1{\mathrm{\Pi }}_{sz}P\left\{y\in \frac{Y_z}{h_s}\right\},$$

where ${\mathrm{\Pi }}_{sz}$ is the loss function, $P\left\{y\in Y_z/h_s\right\}$ is the conditional probability of falling within the domain $Y_z$, if the state $h_s$ is the case, $s\ne z,s,z$ are the state indexes (symbol “/” denotes conditional probability). By averaging the conditional risk function over all the states $h_s$ we obtain average risk

$$R={\displaystyle \sum }_{s=0}^1{p}_s{R}_s,$$

where $p_s$ is the prior probability of the state $h_s$.

For a simple loss function ${\mathrm{\Pi }}_{sz}=\left\{\begin{array}{c}1, s\ne z,\\ 0,s=z,\end{array}\right.$, using posterior probabilities $P\left\{h_s/y\right\}, s=0, 1$, we obtain the posterior risk

$$R={{\displaystyle \sum }}_{s\ne z}P\left\{h_s/yϵY_z\right\}.$$

Minimizing the risk $R$, we obtain the best estimates when there is no a priori knowledge on a useful signal and the presence of a high noise level. Figure 3 shows the results of the application of operations (8) to geomagnetic field data of the Paratunka observatory during the magnetic storm on 12 May 2021. The results show that the application of adaptive thresholds $T_{i,s}$ (Figure 3b) allows for the detection of weak (threshold $T_{1,s}$) and strong (threshold $T_{2,s}$) geomagnetic disturbances before and during a magnetic storm. As a comparison, Figure 3c illustrates the results of the application of the threshold $T=\sigma \sqrt{2\ln N}$ [45]. It was proved in the paper [45] that in white noise case with dispersion ${\sigma }^2$, the threshold $T=\sigma \sqrt{2\ln N}$ allows us to obtain the estimates close to optimal ones. The results show that the application of the threshold $T$ (Figure 3c) does not allow one to localize the features and to suppress noise in a signal completely. The results confirm the efficiency of the suggested method and agree with the results of the paper [46], in which it was shown that the application of adaptive thresholds allows one to obtain more accurate estimates in the case of data with features of a different structure and a high level of noise. The results also confirm the author’s statement [42] that when signal energy is relatively low noise energy, optimal estimates give too small a threshold and its application does not make it possible to obtain satisfactory results.

3. Calculation of Results and Discussion

In the calculations, we used the magnetic data with the 0.5-s rate (raw and pre-processed) obtained at the Paratunka observatory since 2009 and one-minute data calculated according to INTERMAGNET standards.

The choice of the observatory was determined by the presence of representative data sampling and the detailed description of measurement conditions, such as noise and its sources in raw data, and the magnetologists’ qualification level required for the correct processing of the results of the observations. Information on noise was used to estimate the possibility to apply the method in the mode close to real time. The method application results are shown in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10.

3.1. Approximation of Geomagnetic Field Quiet Variations

Figure 4 illustrates the Sq-curves for March and September 2021. They were calculated by a traditional technique (black curve at Figure 4) and based on smoothed components of the 6-level decomposition (green curve). The optimum decomposition level $m^{*}=6$ is determined with the algorithm described above. The results show that Sq-curves obtained by the proposed method allow us to get a better approximation to the quiet variation, compared with the traditional approach, and confirm the method efficiency.

3.2. Detection of Short-Period Geomagnetic Disturbances, Second-Resolution Data Processing

Figure 5, Figure 6 and Figure 7 show the processing results for the raw second data and the data used during the preparation of INTERMAGNET preliminary minute values. The data were obtained on 5 December 2017, by a fluxgate magnetometer FGE-DTU [47], installed at the Paratunka observatory within the scientific cooperation between IKIR FEB RAS and the Deutsches GeoForschungsZentrum GFZ (Helmholtz-Zentrum Potsdam). The magnetometer sensitivity is up to 0.01 nT and the measurement frequency is 2 Hz. A suspended system of sensors is used in the magnetometer. It automatically aligns the Z-axis to a vertical position; however, it causes sensitivity to mechanic impacts on the magnetometer basis, for example, during seismic wave passing or severe shocks near an installation place. There is some noise (oscillations) in the interval under consideration in the records. It is caused by the earthquake that occurred at 09:06:04 UT at a depth of 40 km and about 170 km away from the observatory. The earthquake class was Ks = 10.6 (Figure 5, Figure 6 and Figure 7, red oval). Noise effects from earthquakes can reach one-hundreds of nT depending on distance, depth, and magnitude. A known and quite frequent problem of the observatories carrying out different kinds of measurements is interference from radio devices, for example, from ionosondes and radars located near magnetometers. At the Paratunka observatory, short packages of noise are observed every 15 minutes in records of FGE-DTU magnetometer. They occur during the ionosphere vertical sounding by automatic ionosonde Parus (Figure 5, Figure 6 and Figure 7, green box). One more common issue is the mutual influence of closely located magnetometers. For example, magnetometers with coil systems, such as dIdD GSM-19FD or POS-4, create additional magnetic fields around a sensor during current overlapping in coils. Additional fields also occur around the connecting cables and proton magnetometer sensors during the polarization process. An example of the latest effect is illustrated in Figure 5, Figure 6 and Figure 7, panel (a) Overhauser magnetometer GSM-90 generates noise in FGE-DTU measurements in the form of a meander with the period of 5 s and amplitude of 0.5 nT and more. A detailed description of the noise mentioned above is given in the paper [48].

Standard processing of FGE-DTU records by a magnetologist at the Paratunka observatory includes some operations. The effect from GSM-90 with 5-s cycle is removed by selecting the FGE data out of the time interval when GSM-90 polarization occurs. In this case, the uniform set of values with the frequency of 2 Hz is uniformly gaped with a frequency decrease to 0.4–0.5 Hz. The noise from the ionosonde is regular and occurs at known times, thus, it is removed automatically when a magnetologist has determined sounding parameters and noise duration during sounding. Irregular noise, for example, from earthquakes, is visually identified by a magnetologist from the raw curves. Noised intervals are marked and are not used in further processing [49].

Our algorithms applied in the analysis need the series with uniformly distributed data without gaps. Thus, the gaps were preliminarily filled on a uniform greed by the nearest neighbor method.

Figure 5 shows magnetic field D-component variations (Figure 5a presents raw data, Figure 5d presents data processed by a magnetologist) and the results of method application (Figure 5b,c,e,f; (8), (9) operations were used). The analysis shows evident noise generated by an earthquake and the ionosonde (Figure 5a, red oval and green box). Comparison of the method results obtained for the raw data (Figure 5b,c) and preliminarily processed data (Figure 5e,f), shows a significant impact of noise on the processing results at the times of its occurrences. Evidently, the anomalies occurred as a result of high noise amplitude. We should note that noise manifestation in the method results is observed only near the signal containing noise, compared, for example, to Fourier methods. The result is provided by wavelet finiteness [29,30].

Figure 6 shows the vertical Z-component variations (Figure 6a presents initial data, Figure 6d presents pre-processed data) and the results of method application (Figure 6b,c,e,f). We can see on the graphs that short (duration is several seconds) noise impact from the ionosonde (Figure 6a, green box) was removed after preliminary processing (Figure 6e,f). Noise from the earthquake at 09:07 UT does not have a visible reflection in the time series (Figure 6a, red oval). However, after the preliminary processing (Figure 6d), a false anomaly appeared due to the removing of part of the data (Figure 6e,f), the correction for the nonorthogonality of H-, D-, and Z-sensors is taken into account in Z-variations. Therefore, data removal in variations H and D causes the gap in the Z record. We should note that the problem of correct and effective filling of the gaps in geomagnetic data is topical, does not have indisputably recommended methods at the present time, and is often solved individually for each case.

Analysis of horizontal component H variations (Figure 7) shows that the considered noise has small amplitude and does not significantly impact method results. Regular interference in the shape of a meander also does not affect the processing results. That allows us to apply the suggested method to analyze geomagnetic data (H-components) in online mode and to obtain results with admissible accuracy in case of short low-amplitude noise. Sudden changes caused by external impact (ionosonde or earthquake effects) of large amplitude should be identified and removed by a magnetologist during the preliminary analysis and processing of raw data, or they should be separately marked in processing results. Estimates showed that if noise amplitude is more than 1–3 nT (depending on regular oscillations caused, for example, by GSM-90 interference), they should be removed to obtain reliable results. Filling long gaps in geomagnetic data, especially when choosing inadequate filling methods, causes significant distortions of natural variation changes and the appearance of false anomalies in results.

Figure 8 illustrates the results of method application (operations (8) and (9)) for the beginning of a magnetic storm on 21 April 2017. In operation (8), coefficient values $U_1=1.5$ (threshold $T_{1,s}$) and $U_2=2.5$ (threshold $T_{2,s}$) were used. The moving time window length corresponded to one hour and was $l=7200$ counts for the initial data and $l=2880$ for the preliminarily processed data. The Paratunka observatory data processing results are illustrated in Figure 8d,e (initial data) and Figure 8g,h (data preliminarily processed by a magnetologist). For comparison, Figure 8j,k shows the results of the processing of a raw dataset of the Magadan observatory of IKIR FEB RAS (IAGA code is MGD, 60.05 N, 150.72 E). These data were not preliminarily processed since noise amplitude did not exceed 0.5 nT. To estimate the state in the near-Earth space, the upper part of Figure 8 represents solar wind speed (SWS) (Figure 8a) and interplanetary magnetic field vertical component IMF Bz (Figure 8b).

According to the resource [50] data, SWS increased from 360 to 500 km/s in the middle of the day on 20 April, IMF Bz fluctuations grew up to Bz = ±12 nT. The magnetic storm began on 21 April at 15:19 UT (red vertical line). The analysis of method results (Figure 8d,e,j,h,j,k) shows short-period geomagnetic disturbances at the Paratunka and Magadan observatories before the event. Disturbances increased as the storm beginning was approaching. We can also note the short-time insignificant increase of geomagnetic activity from 13:45 to 14:00 UT at the times of SWS growth (Figure 8a) and insignificant increase of southward IMF Bz fluctuations (Figure 8b). At the time of the storm beginning, disturbance intensity increased significantly. A comparison of the results for the data from different observatories shows clear common dynamics of the process. Still, at the Magadan observatory, the disturbances have the highest amplitude associated with the site location near auroral latitudes. A comparison of the results for initial (Figure 8d,e) and preliminary processed data (Figure 8g,h) shows that noise does not affect the method results significantly and confirms its efficiency.

3.3. Detection of Short-Period Geomagnetic Disturbances, Minute-Resolution Data Processing

Minute-resolution data were obtained using the method defined by INTERMAGNET standards ([51], Section 2.4). Filtering, with Gaussian weight coefficients, was applied to the second data on the interval of −45 … +45 s. Minute value was centered at the beginning of a minute. When some second values were missing, weight renormalization was carried out.

Figure 9 and Figure 10 show the results of the method applied for a moderate magnetic storm on 12 May 2021. To estimate the near-Earth space state, the upper panels of Figure 9 show SWS (Figure 9a), IMF Bz (Figure 9b), and geomagnetic activity Dst-index (Figure 9c). According to the resource [50] SWS was about 300 km/s on 11 May, IMF Bz fluctuations were within ±2 nT. The application of operations (8) and (9) to the PET minute data are illustrated in Figure 9e,f. In operation (8) threshold coefficient $U_2=2.5$ was used (positive and negative disturbances exceeding the threshold $T_{2,s}$ are illustrated in Figure 9 by pink and by blue, respectively). Based on the ratio (8), application of the coefficient $U_2=2.5$ allows one on each scale $s$ to detect wavelet coefficient fluctuations exceeding in amplitude the 2.5 RMSD from their characteristic values ($T_{2,s}=2.5\times S{t}_s$) within the time window under analysis (window length $l=720$ counts, corresponding to 12 hours, was used). To detect weaker geomagnetic disturbances, threshold coefficient $U_2=1.5$ was used (positive and negative disturbances exceeding the threshold $T_{1,s}$ are illustrated in Figure 9 by yellow and light blue, respectively). We should note that the times of weak geomagnetic disturbance occurrences before the event correlate with the corresponding small IMF Bz fluctuations times. Moreover, a clear increase in disturbance intensity is observed as the magnetic storm beginning is approaching. A similar character of pre-storm short-period geomagnetic disturbances, according to the meridional station network and auroral zone stations, was mentioned in the papers [31,32]. Strong magnetic storms on 7 January and 17 March 2015 [32] and weak magnetic storms on 9 July and 27 September 2017 were under analysis [31].

Based on space weather data [50], an inhomogeneous accelerated flux from a coronal hole and coronal mass ejection (CME from 9 May) arrived at 05:49 UT on 12 May. During that period, SWS suddenly increased (Figure 9a), and IMF BZ fluctuations increased to Bz = ±20 nT (Figure 9b). The time of the storm’s sudden commencement is evident on the magnetic field H-component graph (Figure 9d, red vertical line). The results of data processing (Figure 9e,f) show that geomagnetic disturbances exceeding the threshold $T_{2,s}$ occurred at the Paratunka observatory from 00:00 UT (12:00 LT) to the time of storm beginning on 12 May. Several hours before the storm, during SWS sudden increase (Figure 9a) and IMF Bz southward turn (Figure 9b), disturbance intensity increased significantly and reached $I_k{}^{2\pm }\left(k\right)=4000$ (Figure 9f). That is likely to be associated with the beginning of compression region formation.

Strong geomagnetic disturbances were observed during the IMF Bz southward turns and fluctuation intensification during the storm initial phase. IMF Bz northward direction accompanies a geomagnetic activity decrease. During the IMF Bz long southward turn, the storm main phase occurred. It was characterized by a sudden decrease in Dst-index (Figure 9c) and significant increase in H-component fluctuations (Figure 9d). We should note that the Paratunka observatory was in the night sector (00:00–02:00 LT) during that period. Based on the processing data, the intensity reached the values $I_k{}^{2\pm }\left(k\right)=12,000$ near the H-component minimum values. That is likely to be associated with ring current effects and intensification of the magnetosphere tail currents. The results agree well with those represented in the papers [52,53].

Figure 10 shows the processing results for the Paratunka observatory data for the same period, but the processing was sequential up to the storm beginning. The obtained results of processing (Figure 10, operations (8), (9) applied) agree with the results illustrated in Figure 9. That confirms the method efficiency and shows the possibilities of its application for fast detection of short-period geomagnetic disturbances preceding a magnetic storm beginning.

4. Conclusions

Investigation results confirmed the proposed method’s efficiency for detecting short-period geomagnetic disturbances of differing intensity and duration. Adaptivity of the wavelet model and high time resolution of the suggested approach allow us to investigate sudden small-amplitude variations of the geomagnetic field, which are of interest in space weather. The method makes it possible to process the data of different resolution and, compared to neural networks, does not require periodic rearrangement of algorithm parameters.

Comparison of the proposed wavelet model with the Sq-curve method showed its advantage. The model allows one to obtain a better approximation to geomagnetic field quiet variation than the traditional approach. The method efficiency was also confirmed on the basis of the comparison of its results with the results of application of the optimal threshold suggested in the paper [45].

Comparing the processing results of raw geomagnetic data with the data preliminary processed by a magnetologist confirmed the possibility of applying the method in online mode. In the case of short low-amplitude noise (amplitude is 1–3 nT), the method makes it possible to obtain results with admissible accuracy and can be used to perform data online analysis. The result is important as long as a great and constantly growing volume of recorded geomagnetic data cannot be analyzed manually within a reasonable amount of time. Sudden high amplitude changes caused by external impacts (for example, ionosonde or earthquake effect) can generate false anomalies. In that case, preliminary processing of the data by a magnetologist is required to improve the analysis accuracy. That indicates the necessity to develop automatic detection methods and compensation of noise in geophysical data. One of the possible solutions to the problem was considered in the papers [6,54].

An analysis of the data during disturbed periods confirmed the investigation results [31,32] and showed the possibility of short-period geomagnetic disturbances preceding a geomagnetic storm beginning. Weak short-period increases in geomagnetic activity before storms were described in a number of papers [3,4,5,55]. During the events under consideration, geomagnetic disturbance intensity increased during SWS growth and intensification of southward IMF Bz. As the storm beginning approached, disturbances were intensified. The observed correlation of the detected geomagnetic disturbances with interplanetary environment parameters indicates their external nature and relation with an oncoming magnetic storm.

Summarizing the work results shows the complicated dynamics of the processes in the near-Earth space during increased solar activity and magnetic storms. Investigation of the processes requires the development of methods for data analysis and an extended observation network. The results confirmed the importance of ground observations, which, together with space data, are a valuable and informative tool for space weather monitoring. The investigation also showed the essential role of site observations, which significantly contribute to the investigation of near-Earth space processes and space weather forecast.