Method of Wavelet-Decomposition to Research Cosmic Ray Variations: Application in Space Weather

Abstract

Since their discovery, cosmic rays have been an integral part of the development of fundamental physics, from the discovery of radiation coming to the Earth from outer space and the identification of high-energy particles in it, as well as new fundamental symmetries in the laws of nature, to the knowledge of residual matter and magnetic fields in interstellar space. Cosmic rays are used in a number of fundamental and applied research in solar-terrestrial physics and are important in the research of the near-Earth space processes. Cosmic ray variations observed on the Earth’s surface are an integral result of various solar, heliospheric, magnetospheric and atmospheric phenomena. The most significant changes in cosmic ray parameters are caused by coronal mass ejections and subsequent changes in the parameters of the interplanetary magnetic field and solar wind. Therefore, the study of cosmic rays makes it possible to obtain valuable information about the processes in the near-Earth space and in the Earth’s magnetosphere during disturbed periods. This article proposes a method for analyzing cosmic ray variations. It is based on the use of wavelet data decomposition operations and their combination with threshold functions. By using adaptive thresholds, the operations for detecting anomalous changes in data and for suppressing the noise were developed. Anomalies in cosmic rays can cause radiation hazard for astronauts, radio communication failures, as well as malfunctions in satellites, leading to the loss of orientation and destruction. Therefore, the task of timely diagnostics of anomalies is urgent. The paper describes the algorithms for the implementation of the method and shows their application in the space weather problem. We used data from the network of ground stations of neutron monitors. The efficiency of the method for detecting abnormal changes of different amplitudes and durations is shown. Application of the method made it possible to detect clearly and to evaluate Forbush effects in cosmic rays, which precede the onset of magnetic storms of various nature and strength.

1. Introduction

Cosmic rays are an integral part of the development of fundamental physics. From the discovery of radiation coming to the Earth from outer space and the identification of high-energy particles in it, as well as new fundamental symmetries in the laws of nature, to the knowledge of residual matter and magnetic fields in interstellar space [1]. Cosmic rays are used in a number of fundamental and applied research in solar-terrestrial physics and are important in the research of the near-Earth space processes. During anomalous processes on the Sun and magnetic storms occurring in the near-Earth space, the possibility of conducting a real time analysis of geophysical parameters variations with acceptable accuracy becomes especially relevant. The difficulty of realizing this possibility is due to the complex, nonlinear structure of the registered data, a wide range of anomaly patterns, as well as a high level of noise of various nature [1,2,3]. During strong magnetic storms, strong changes occur in radiation belts. Deformation of the magnetosphere and, as a consequence, compression of the radiation belt and displacement of its maximum extent are observed [4,5]. These deformations are characterized by an increase in the number of low-energy electrons and a decrease in high-energy particles, and then an increase in the latter [4]. These changes occur, as a rule, within a day. During the reverse phase of a magnetic storm, the intensity of cosmic radiation increases, sometimes exceeding the intensity before the storm. Such events confirm the existence of particle acceleration mechanisms during magnetic disturbances [5,6]. Within a few weeks, everything is usually recovered, but sometimes irreversible changes occur, and a large number of fields inhomogeneities arise, causing particle diffusion deep into the magnetosphere [7]. In their turn, rapid variations in the outer part of the belts indicate the existence of efficient high-speed mechanisms for the replenishment of the outer belt with electrons [4]. Therefore, cosmic rays (CR), which are the object of this research, are one of the significant indicators in space weather.

Correlated observations of cosmic rays on a global scale remain insufficiently studied [8,9,10,11,12,13,14,15]. Cosmic ray variations are divided into periodic and non-periodic. Periodic ones include 22-year, 11-year, 27-day and solar-daily variations. Non-periodic variations (anomalous, sporadic) include Forbush effects (sharp sudden change in the intensity of cosmic rays) and GLE-events (ground level enhancement of solar cosmic rays) [16]. Anomalous variations in cosmic rays can cause radiation hazard for astronauts, failures in radio communications, as well as malfunctions in satellite operation leading to loss of orientation and destruction. Therefore, the problem of timely diagnosis of such events is topical and is being solved by a large number of scientific groups [8,9,10,11,12,16,17,18,19,20,21].

A preliminary increase or decrease in CRI can be observed before a magnetic storm. That makes it possible to use it as a precursor (for example, [20,21]). However, CR anomalies observed on the Earth may be disassociated with a magnetic storm [22,23]. Therefore, to improve the accuracy of predicting magnetic storms, it is important to identify anomalies in CR. This requires research involving various possible methods. The investigation described in the paper [15] suggests that all CME properties show some correlation with CRI, but the question of determining the parameter that has the highest correlation with CRI, is still open. In [11] the authors studied the peculiarities of the cosmic ray intensity behavior during a strongly disturbed period of 4–10 September 2017. Analysis of the correlation between the CRI and Dst variability during that strong storm showed the presence of a Dst delay of 3–4 h. This indicates the importance of taking into account the CR variability measure along with the IMF and solar wind characteristics in space weather. It is important to note that the correlation of the initial data of the CR intensity with Dst obtained by the authors [11] had the maximum value with a delay close to zero. Therefore, the initial CR data in space weather are ineffective and it is necessary to use more accurate and sensitive methods of analysis.

At present, cosmic ray analysis is performed using different methods [8,9,10,24]. The use of the global survey method (GSM) [25] for calculating spherical harmonics is not effective for all anomalous events in cosmic rays (for example, beginning of the proton increase) [26]. The second disadvantage of GSM is the dependence on the number of cosmic ray observation stations that are definitely located around the globe and that are not always in good working order. The third disadvantage of this method is the complexity and duration of the calculations, which do not allow one to receive a prompt result. The applied averaging methods [27] make it possible to study the characteristic changes in GSM, but they are insensitive to low-amplitude anomalies that appear on the eve of strong solar flares or Forbush effects. Another disadvantage of the averaging methods is the risk of information distortion, that can result both in an undetected anomaly and a “false” alarm. The research conducted by authors [28,29,30] has shown that the use of machine learning methods can effectively detect multiscale anomalies in cosmic rays of a narrow spectrum, but when different-scale anomalies occur, the efficiency decreases. Moreover, these methods are also sensitive to signal changes associated with the solar activity level. In a number of works [11,12,31,32,33], it was shown that the wavelet transform apparatus is the most effective for describing non-stationary processes containing multiscale characteristics. Wavelet transform allows one to detect complex structures in data and to discover singularities and transient phenomena [11,12,31]. While studying the cosmic ray quasiperiodic variations, the authors [12], compared the wavelet transform with the classical FT method (the spectral power density was estimated). The obtained results confirmed the effectiveness of studying the GSM time evolution and revealing quasiperiodic features using the wavelet transform. For the detection of low-amplitude periodicities, the authors [12] used an F-filter [34], which allows one to assess the variability of the process within a given time interval (using a time window) and to detect the “hidden” patterns in data more clearly.

This article proposes a method for analyzing cosmic ray variations. It is based on a heuristic approach and is a combination of the method that we proposed for constructing an adaptive nonlinear approximating scheme in the basis of orthogonal functions [35,36] with an anomaly detection algorithm [33] using the wavelet spectrum construction [35,37,38]. Using the data of neutron monitors, the method allows us to identify anomalous changes in the rate of galactic cosmic ray arrival on the Earth. This approach is close to the approach used by the authors in the paper [12] and allows detecting multiscale features of the CR intensity behavior. Following the results of the work [28], to enhance the effect of detecting small anomalies, the paper proposes to use nonlinear adaptive threshold estimates before applying continuous wavelet transform. By the example of the basis sets of wavelet packets, the high sensitivity of the method and its effectiveness for detecting small anomalies in cosmic rays during weak and moderate disturbances in the near-Earth space were shown. The problem of detecting an anomaly of weak intensity is the most difficult and requires high sensitivity of methods, as well as the possibility to make real time calculations. As noted above, the significance of the problem of detecting such anomalies is determined by their occurrence before strong disturbances in the magnetosphere and serving as their predictors [9,16]. Since the anomalous variations in cosmic rays at different observation stations have different structures and intensities or may be absent, it is important, as was noted by the authors [8], to analyze the data from a station network.

The paper presents the results of the method using the data from ground-based neutron monitors Inuvik (INV, Coord: 68.36, −133.72), South Pole (SOPO, Coord: −90, 0) and Thule (THUL, Coord: 76.5, −68.7) [39]. Structures in the data were detected. They have a complex spectrum, are correlated in space and time, and are recorded sequentially with a time delay from 2 to 6 h at Inuvik and Thule stations. The results showed the possibility of using the method for a detailed study of cosmic ray flux dynamics using the data of neutron monitors and the timely detection of anomalies of different structure and intensity.

2. Materials and Methods

2.1. Identification of Informative Signal Structures and Noise Reduction

Following the results of the work [28], the signal $f\in H$ ($H$ is the Hilbert space) is represented by the vectors adaptively selected from the orthonormal basis $\mathfrak{B}={\left\{g_m\right\}}_{mϵ\mathbb{N}}$ ($\mathbb{N}$ are natural numbers, including 0) of the space $H$:

$$f_I={\sum }_{m\in I}\langle f,g_m\rangle g_m,$$

(1)

where $I$ is the set of indices.

The approximation error is estimated as $ϵ\left[I\right]=\Vert f-f_I\Vert {}^2={\sum }_{m\notin I}{\left|\langle f,g_m\rangle \right|}^2$.

In order to minimize the error $ϵ\left[I\right]$, we choose $I$ in such a way that the vectors $g_m$с with indices from $I$ have the largest moduli of the scalar product $\left|\langle f,g_m\rangle \right|$. In this case, representation (1) can be obtained by applying the threshold function [35]

$$T\left(x\right)=\left\{\begin{array}{c}x,if\left|x\right|\ge T,\\ 0,if\left|x\right|<T.\end{array}\right.$$

Then we have

$$f_I={\displaystyle \sum _{mϵ\mathbb{N}}}T\left(\langle f,g_m\rangle \right)g_m.$$

(2)

For discrete and noisy signal $f\left[m\right]$:

$$X\left[m\right]=f\left[m\right]+V\left[m\right],$$

(3)

where $X\left[m\right]$ is the recorded discrete data, $f\left[m\right]$ is an informative signal, $V\left[m\right]$ is noise, from (2) we obtain the estimate

$$\tilde{F}=DX={\displaystyle \sum _{mϵ\mathbb{N}}}T\left(\langle X,g_m\rangle \right)g_m.$$

(4)

The risk of estimating $\tilde{F}$ is $r\left(D,f\right)=E\left\{\Vert \tilde{F}-f{\Vert }^2\right\}$, where $E$ is the mathematical expectation.

To minimize the risk $r\left(D,f\right)$ we can use the minimax approach [40], and the problem is to determine an operator D that minimizes the maximum risk (minimax risk) [40]: $r_O\left(ϴ\right)=in{f}_{DϵO}su{p}_{fϵϴ}E\left\{\Vert \tilde{F}-f{\Vert }^2\right\}$, $O$ is a set of operators performing mapping (4).

To minimize the risk $r_O\left(ϴ\right)$ the threshold $T$ is chosen so that there is a high probability that it is greater than the maximum level of noise figures $\left|V_{\mathfrak{B}}\left[m\right]\right|,V_{\mathfrak{B}}\left[m\right]=\langle V,g_m\rangle$ (see (3)) [35]. We also take into account that the risk of estimate (4) is associated with the approximation error f in basis $\mathfrak{B}$.

Using the bases of wavelet packets [35,38] ${\mathfrak{B}}_j^p={\{{\Psi }_j^p(2^{j}t-m)\}}_{m\in \mathbb{N}}$, where ${\mathfrak{B}}_j^p$ is the basis of the space $W_j^p$, on the basis of mapping (4), we obtain a tree decomposition of wavelet packets. In the decomposition, the space $W_j^p$ is divided into orthogonal subspaces

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

(5)

Following the work [28], in each node of the tree we use an adaptive threshold $T_j$ (the formula for determining the threshold is given in the algorithm below) that with a high probability exceeds the maximum level of noise factors $\left|V_{{\mathfrak{B}}_j^p}\left[m\right]\right|$, ensures its suppression and maintains the coefficients in the vicinity of structural features of the signal [35]. The algorithm for obtaining the estimate (4) (ACAS algorithm [28]) in this case includes the following operations:

(1) We decompose the signal $X$ into wavelet packets (see (5)):

$$W_j^0:W_j^0={\oplus }_{i=0}^I{W}_{j_i}^{p_i},{\{{\Psi }_{j_i}^{p_i}(2^{j_i}t-m)\}}_{m\in \mathbb{N}}$$

is the basis of the space $W_{j_i}^{p_i}$;

(2) Based on the estimate of normalized energies, we determine the tree branches corresponding to the structural components of the signal: the basis $B_{j_i}^{p_i}$ of the space $W_{j_i}^{p_i}$ is the basis:

$$B_{j_i}^{p_i}=\left\{\begin{array}{c}{\{{\Psi }_{j_i}^{p_i}(2^{j_i}t-m)\}}_{m\in \mathbb{N}},if{\displaystyle \sum _{m\in I^{p_i}}}{\left|\langle X,{\Psi }_{j_i,m}^{p_i}\rangle \right|}^2\ge {\displaystyle \sum _{m\in I^{2p_i}}}{\left|\langle X,{\Psi }_{j_i+1,m}^{2p_i}\rangle \right|}^2+{\displaystyle \sum _{m\in I^{2p_i+1}}}{\left|\langle X,{\Psi }_{j_i+1,m}^{2p_i+1}\rangle \right|}^2\\ {\left\{{\Psi }_{j_i+1}^{2p_i}\right\}}_{m\in \mathbb{N}}\cup {\left\{{\Psi }_{j_i+1}^{2p_i+1}\right\}}_{m\in \mathbb{N}},if{\displaystyle \sum _{m\in I^{p_i}}}{\left|\langle X,{\Psi }_{j_i,m}^{p_i}\rangle \right|}^2<{\displaystyle \sum _{m\in I^{2p_i}}}{\left|\langle X,{\Psi }_{j_i+1,m}^{2p_i}\rangle \right|}^2+{\displaystyle \sum _{m\in I^{2p_i+1}}}{\left|\langle X,{\Psi }_{j_i+1,m}^{2p_i+1}\rangle \right|}^2\end{array}\right.,$$

(6)

where the set of indices $I^l,l=p_i,2p_i,2p_i+1$ are defined as follows: index $m\in I^l$, if $\left|\langle X,{\Psi }_{j_i,m}^l\rangle \right|\ge T_{j_i}$, threshold $T_{j_i}=K\times {\sigma }_{j_i}^l$, ${\sigma }_{j_i}^l=\sqrt{\frac{1}{L}{\sum }_{m=1}^L{(\langle X,{\Psi }_{j_i,m}^l\rangle -\langle \overline{X,{\Psi }_{j_i,m}^l}\rangle )}^2}$, where the threshold coefficient K is determined by estimating the posterior risk, $\overline{\langle X,{\Psi }_{j_i,m}^l\rangle }$ is the mean of the set ${\left\{\left|\langle X,{\Psi }_{j_i,m}^l\rangle \right|\right\}}_{0\le m<L}$, $L$ is the number of elements.

The nodes of the tree ($j_i$, $p_i$) selected on the basis of Operation (6) define informative signal components. Figure 1 shows the trees constructed according to the data of stations monitors Inuvik, Thule and South Pole. They indicate the presence of differences in the data structure at different latitudes and a wide range of -patterns of the forming features. Following the results of the work [28], we used the Coiflet 2 [38] basis and applied the threshold coefficient K = 2.5 (Operation (6)) when processing the data. The result confirms the complex structure of cosmic ray data and the presence of high noise levels. Comparison of the initial data (Figure 1b,e,h) with the data obtained after the processing (Figure 1c,f,i) shows the effectiveness of the proposed algorithm.

2.2. Anomaly Detection and Estimate of Their Intensity

(1) Based on the BAS algorithm, we obtain a signal representation:

$${\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}=D_{{\mathfrak{B}}_{j_i}^{p_i}}X={\displaystyle \sum _{j_i}\sum _m}{T}_{j_i}\left(\langle X,{\Psi }_{j_i,m}^{p_i}\rangle \right){\Psi }_{j_i,m}^{p_i},$$

where ${\mathfrak{B}}_{j_i}^{p_i}={\{{\Psi }_{j_i}^{p_i}(2^{j_i}t-m)\}}_{m\in \mathbb{N}}$ are wavelet packet bases; tree nodes ($j_i$, $p_i$), defining informative signal components.

Using the continuous wavelet transform, we obtain the wavelet spectrum of the resulting map [38]:

$${\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}W{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\left(u,s\right)={\int }_{-\infty }^{+\infty }{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\frac{1}{\sqrt{s}}\Psi \left(\frac{t-u}{s}\right)dt,$$

where $\Psi$ is a wavelet, $s$ is a scale, $u$ is a time shift, $s,u\in ℝ$ ($ℝ$ are real numbers), $s\ne 0$.

(2) The amplitude of the coefficients $\left|W{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\left(u,s\right)\right|$ characterizes the amplitude of the local feature of the function on the scale $s$ in the vicinity of the point $t=u$ [41], therefore, the amplitude increase indicates the anomaly occurrence in the vicinity of this point. To detect anomalies on the $s$ scale, we use the thresholds $T_s^l$, which, taking into account the change in the dynamics of the process, are calculated in a sliding time window:

$$P_{T_s}\left[W{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\left(u,s\right)\right]=\left\{\begin{array}{c}W{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\left(u,s\right),if\left|W{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\left(u,s\right)\right|\ge T_s^l,\\ 0,if\left|W{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\left(u,s\right)\right|<T_s^l,\end{array}\right.$$

(7)

where $T_s^l$ = q × ${\sigma }_s^l$, ${\sigma }_s^l$ is the standard deviation (SD) of the coefficients, calculated in a sliding window of length $l$, $\mathrm{q}$ is the threshold coefficient.

(3) For each moment of time $t=u$ we estimate the intensity of the anomaly

$$E_u={\displaystyle \sum _s}{P}_{T_s}\left[W{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\left(u,s\right)\right],$$

(8)

which in the case of an anomalous increase (relative to the characteristic level) will have a positive value and, in the case of an anomalous decrease, will have a negative value.

Figure 2 shows an example of applying the proposed method to the neutron monitor data of the Thule station for 26–28 May 2019. Figure 2b shows the result of the application of the ACAS. When performing operation (7), the thresholds $T_s^l=$2$\times {\sigma }_s^l$ (Figure 2c,e,g,i) and $T_s^l=$2.5$\times {\sigma }_s^l$ (Figure 2d,f,h,j), were used. The standard deviation ${\sigma }_s^l$ was estimated in a sliding time window of the length $l=1440$ that corresponds to a day. The estimation of the intensity of the anomalies $E_u$ (Operation (8)) was performed in two ways:

$$E_{u,1}={\sum }_{s=1}^{2600}{P}_{T_s}\left[W{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\left(u,s\right)\right],E_{u,2}={\sum }_{s=1000}^{2600}{P}_{T_s}\left[W{\tilde{F}}_{{\mathfrak{B}}_{j_i}^{p_i}}\left(u,s\right)\right].$$

The positive values of the intensity $E_u$ are shown in Figure 2e,f,i,j in red, negative ones are shown in blue. According to IZMIRAN data [42], in the presented period, Forbush effects were recorded at 22:14:00 on 26 May 2019 (MagnM = 1.1) and at 15:00:00 on 28 May 2019 (MagnM = 0.5). In Figure 2, they are marked with red lines. According to IPG data [43], geomagnetic disturbances at high latitudes were registered at 3:00:00 on 27 May 2019 and at 19:00:00 on 28 May 2019. In Figure 2, they are marked with yellow lines. The results of data processing show that there was a gradual increase in the intensity of cosmic rays at the analyzed station during the first Forbush effect. The increase exceeded the background level 23 h before geomagnetic activity growth. During the disturbances, the cosmic ray intensity decreased and reached the background level. The second anomalous period is characterized by a decrease in cosmic ray intensity that occurred 18 h before the geomagnetic disturbances. Further, a slight increase in cosmic ray intensity was observed. The results show complex dynamics of cosmic rays during the disturbed periods.

Comparison of the processing results when using different thresholds shows that a higher threshold $T_s^l=$2.5$\times {\sigma }_s^l$ (Figure 2d,f,h,j) localizes anomalous feature in the signal more clearly. We should also note that the calculation of the intensity by the $E_{u,2}$ method (scales from 1000 to 2600 are used, Figure 2g–j) allows us to reduce the influence of uninformative small-scale variations associated with daily variation, precipitation, etc. Note that the identified anomalous changes in cosmic rays had small amplitude an, based on the analysis of the initial variations (Figure 2a), it is very difficult to distinguish them. This indicates the high detecting ability of wavelets and the efficiency of the proposed method. The effectiveness of the method also lies in the ability to detect anomalies automatically (without an expert).

The identified anomalies in cosmic ray variations are consistent with the results of the studies presented in the papers [9,10,16,37], and confirm the possibility of their occurrence before and during geomagnetic disturbances.

Figure 3. for comparison, shows the results of detecting Forbush effects in NM data based on the proposed method (Figure 3e–h) and the results of using CWT (Figure 3c,d). Figure 3a,b shows data from neutron monitors at Inuvik and Thule South Pole stations. During the analyzed period, the Forbush effect was recorded at 18:00 UT on 3 May 2019 according to IZMIRAN data [42]. It is marked in Figure 3 with a red vertical line. Moments of increased geomagnetic activity at 04:00 UT on 3 May (min Dst = −17, weak storm at high latitudes) and at 04:00 UT on 4 May (min Dst = −22) are marked in Figure 3 with a yellow vertical line. According to the data of [43], due to the arrival of an inhomogeneous accelerated flow from the western part of the coronal hole, the SWS reached 579 km/s on 2 May. The southern IMF component fluctuations increased to $B_z$ = ±10 nT. During that period the results of the method show a weak anomalous increase in the CR intensity at the Thule station (Figure 3h). At the end of the day on 3 May, an inhomogeneous accelerated flow came from the eastern part of the coronal hole, the SWS on May 04 increased to 526 km/s, the southern IMF component fluctuations increased to $B_z$ = ±8 nT [43]. According to the processing data (Figure 3e–h) at the end of the day on May 3, at the time of the Forbush effect, an anomalous decrease in the CR intensity occurred at Inuvik and Thule stations. The decrease is more pronounced at Thule station (it exceeded the background level 4–6 hours before the beginning of the magnetic storm). Until the end of the period, due to the weakening of the coronal hole influence, the SWS decreased to 320 km/s, the southern IMF component fluctuated from $B_z$ = ±4 nT to $B_z$ = ±1 nT. Comparison of the CWT results (Figure 3c,d) with the results of the proposed method (Figure 3e–h) shows its higher accuracy. The CWT results are consistent with the results of the method. However, due to the influence of noise factors (in particular, they are associated with the daily CR variation observed in the CWT spectrum), it is difficult to detect anomalies in CR based on CWT. Moreover, the advantage of the proposed method lies in the possibility of obtaining a more accurate quantitative measure of CR variability (Figure 3g,h).

Evaluation of the proposed method efficiency is presented in

Table 1. Estimation of the proposed method efficiency.

Year	Number of FE/Number of False FEs (No Storm)	Method Results
2013–2015 (high SA)	285/32	Detected: ~86%
		Not detected: ~14%
		False alarm 1 (no FE): ~13% False alarm 2 (no GS): ~9%
2019 (low SA)	97/15	Detected: ~89%
		Not detected: ~11%
		False alarm 1 (no FE): ~11% False alarm 2 (the presence of a FE, but no GS): ~8%

. Table 1 shows that detection of Forbush effects during high solar activity (SA) is ~86%, during low solar activity it is ~89%. The frequency of false alarms during high SA is ~13% (false alarm 1 in Table 1), during low solar activity it is ~11%. The detection rate of anomalies in the CR not related to GS (False alarm 2 in Table 1) is ~9%. The results confirm the method’s effectiveness. To increase the frequency of detecting anomalies in CR preceding GS, it is necessary to attract a larger set of analyzed stations and expand statistics. In addition, an increase in forecast efficiency is obviously possible with the use of other parameters of the interplanetary medium and phenomena on the Sun, based on a comprehensive analysis of the data. Moreover, additional research and deeper study of the processes in the near-Earth space during disturbed periods are required. The variability of cosmic rays and its relationship with solar processes have not been sufficiently studied yet.

3. Results

During the first analyzed period, 24–29 May 2019 (Figure 4), two weak magnetic storms occurred at high latitudes [43]. Figure 4d,e,f shows data from neutron monitors at Inuvik, Thule and South Pole stations. Figure 4g–l shows the results of data processing from the neutron monitors at Inuvik, Thule and South Pole stations. To analyze the state of the near-Earth space, the figure also shows the data on the solar wind speed (SWS) (Figure 4a), values of the interplanetary magnetic field (IMF) $B_z$ component (Figure 4b) and the data on the geomagnetic activity Dst-index (Figure 4c). On the eve of the first event, the near-earth space was calm. Such a state allows us to obtain results with the smallest error, when processing the data in a sliding time window (a daily time window of 1440 counts was used). According to the processing of the data (Figure 4g–l) from 24 May to 25 May, the cosmic ray intensity was within the background values. At about 6:00 UT on 26 May, a smooth increase in cosmic rays’ intensity was observed at all analyzed stations (Figure 4g–l). It exceeded the background level and reached a peak at about 15:00 UT at the stations Inuvik and South Pole. According to the resource [43], that situation coincided with the arrival of a nonuniform accelerated flow from a coronal mass ejection (CME occurred around 21:25 UT on 26 May, the proton density was about 27 particles/cm³). By the end of the day on 26 May, SWS increased (Figure 4a) and remained within 350–400 km/s until 28 May. IMF fluctuations increased to $B_z$ = ±9.5 nT (Figure 4b). According to IZMIRAN [42], at 22:14 UT on 6 May, the Forbush effect was recorded. It is marked with a red line in Figure 4. During that period, positive values of the Dst index were observed (Figure 4c). They are characteristic for the initial phase of the magnetic storm [44] that occurred at high latitudes at 3:00 UT on 27 May (4 h after the recorded Forbush effect). The moment of the magnetic storm beginning is marked with a yellow line in Figure 4. The result is consistent with that obtained in [11].

According to the data processed (Figure 4g–l), the cosmic ray intensity was within the background during the event. Further, from the beginning of the day on 28 May, a smooth decrease in cosmic ray intensity (Forbush decrease) was observed during the period of $B_z$ positive values at the Inuvik and Thule stations (Figure 4g,h,j,k). The Forbush decrease reached its lowest values at the moments of $B_z$ turns to the south and, according to [43], coincided with the arrival of a nonuniform accelerated flow from the coronal hole (CIR). Note that the Forbush decrease had a small amplitude and was more expressed at the Tule station. It exceeded the variation background level at 9:00 UT.

In the second half of the day on 28 May, SWS began to increase (Figure 4a), IMF fluctuations increased to $B_z$ = ±8 nT (Figure 4b). According to IZMIRAN [42], the Forbush effect was observed at 15:00:00 UT. The beginning of the magnetic storm at high latitudes was recorded on 28 May 2019 at 19:00 UT [43]. That coincides with the moment of a sharp southward turn of the IMF $B_z$ component (Figure 4b). During the initial phase of the storm, a slight increase in cosmic rays’ intensity of was observed (Figure 4g–l) at Inuvik and Thule stations. Note, despite of the differences in the data of neutron monitors from different stations, the processing results (Figure 4g–l) show the presence of a clearly expressed general character in the GCR dynamics both before and during the event. That confirms the effectiveness of the method and shows the large-scale nature of the identified anomalous changes. We should also note that the first analyzed magnetic storm was caused by CME, the second one was generated by CIR, and the GCR dynamics before the events was different. The result indicates complex processes in the near-earth space during increased solar activity and magnetic storms, the research of which requires a wide network of observations and development of data analysis methods.

During the next analyzed period, 24–29 August 2021 (Figure 5), two magnetic storms of different strengths occurred. Figure 5d,e,f shows data from neutron monitors at Inuvik, Thule and South Pole stations. At the end of the day on 24 August, a uninform accelerated flux from the coronal hole arrived (CIR), the maximum concentration of protons of which was 27 particles/cm³ [43]. Fluctuations of the southern IMF component increased to $B_z$ = ±7 nT (Figure 5b), SWS increased to 380 km/s (Figure 5a). The moment of the beginning of a weak magnetic storm is marked with a yellow line in Figure 5.

According to the processing results (Figure 5g–l), Forbush decrease occurred at all analyzed stations before the event, during a smooth increase in SWS (Figure 5a). It reached the maximum amplitude at Inuvik station (at 12:00 UT). During the initial phase of the magnetic storm, another Forbush decrease was observed at Inuvik and Thule stations (Figure 5g,h,j,k) at the background of a sharp increase in SWS (Figure 5a) and a prolonged southward turn of $B_z$ (Figure 5b). The Forbush decrease reached its maximum amplitude during the main phase of the magnetic storm and continued until $B_z$ turned northward. The detected effects in cosmic rays are apparently caused by the fast fluxes of the solar wind inside the magnetic cloud (ICME), which shield the GCR with a strong internal magnetic field and lead to the Forbush effect [42].

At the beginning of the day on 27 August, a uninform accelerated flux from a coronal hole and a coronal mass ejection (CME of 23 August) arrived. Its proton concentration was more than 30 particles/cm³ [43]. During that period, an increase in the fluctuations of the southern component up to $B_z$ = ±15 nT (Figure 5b) and an increase in SWS to 545 km/s (Figure 5a) were observed. That led to the occurrence of a strong magnetic storm (minimum Dst = −82) recorded at 02:00:00 UT on 27 August (Figure 5c) [43]. The beginning of the storm is marked with a yellow line in Figure 5.

An anomalous increase in cosmic rays intensity is observed at Inuvik and Thule stations before the event (Figure 5g,h,j,k). It reached the highest amplitude at Thule station 5 hours before the beginning of the magnetic storm. At South Pole station, the GCR flux intensity remained low during that period. At the time of the amplitude increase and a sharp turn of $B_z$ to the south (Figure 5b), a sharp short-term increase in cosmic rays intensity was observed at South Pole station (Figure 5i,l). Note that during the period of the strongest geomagnetic disturbances (the main phase of the storm), short-period fluctuations in cosmic ray intensity were observed at all stations (short-term decreases and increases, Figure 5g–l). These fluctuations indicate interplanetary medium changes and complex nature of the processes.

Similarly to the previous period, the processing results show the presence of a clearly expressed general character in cosmic ray dynamics at all analyzed stations.

4. Conclusions

The performed analysis indicates complex dynamics of the GCR during increased solar activity and magnetic storms. The research of the dynamics requires a wide network of observations and an extensive set of methods and data of geophysical monitoring. The considered events confirmed the possibility of anomalous changes in the GCR during the periods preceding the onset of magnetic storms of different nature and strength. The results confirm the studies of the authors [11,20,21] and show the possibility of using the CRI variability measure as a precursor to a magnetic storm. However, since the CR anomalies may be disassociated with a magnetic storm [22,23], additional studies and a deeper study of the processes in the near-Earth space during disturbed periods are required to improve the forecast accuracy.

The article shows that cosmic ray dynamics at different stations is of a general nature. The Forbush effects can be characterized both by an increase and a decrease in the cosmic ray flux intensity, which reaches its maximum value several hours before the onset of magnetic storms. This confirms the correlation obtained in the paper [11] between the GCR variability and Dst, which has a maximum with a delay of about 3–4 h.

The observed correlation with the changes in the interplanetary medium parameters indicates the reliability of the obtained results. The application of the method makes it possible to detect the Forbush effects that appear at the background of southward turns and an increase in the amplitude of IMF $B_z$ fluctuations, using the data of neutron monitors.

The comparison of the method with CWT also confirmed its effectiveness. The CWT results are consistent with the results of the method, but, in comparison with the CWT, it allows us to detect anomalies in the GCR more clearly and to obtain a more accurate estimate of the CR variability measure.

Evaluation of the method effectiveness for 2013–2015 and 2019 showed that the Forbush effect detection rate of is more than 86%, while the detection of anomalies in CR not related to GS (false alarm) is about 9%. This indicates the need to develop methods for data analysis and space weather forecasting. An integrated approach, using different parameters of the interplanetary medium and phenomena on the Sun, is required.

Thus, the work results confirmed the effectiveness of the method for analyzing cosmic ray data and detecting anomalous changes of different intensity and duration. Using the considered events as an example, the possibility of application of the method for detecting low-amplitude Forbush effects that can precede the onset of magnetic storms and serve as their predictors, has been shown. Numerical implementation allows us to apply the method in operational analysis when making space weather forecasting that determines the applied significance of the research.