**Trapped Proton Fluxes Estimation Inside the South Atlantic Anomaly Using the NASA AE9/AP9/SPM Radiation Models along the China Seismo-Electromagnetic Satellite Orbit**

**Matteo Martucci 1,2,***∗***, Roberta Sparvoli 1,2, Simona Bartocci 1, Roberto Battiston 3,4, William Jerome Burger 4,5, Donatella Campana 6, Luca Carfora 1,2, Guido Castellini 7, Livio Conti 1,8, Andrea Contin 9,10, Cinzia De Donato 1, Cristian De Santis 1, Francesco Maria Follega 3,4, Roberto Iuppa 3,4, Ignazio Lazzizzera 3,4, Nadir Marcelli 1,2, Giuseppe Masciantonio 1, Matteo Mergé 1,†, Alberto Oliva 10, Giuseppe Osteria 6, Francesco Palma 1,†, Federico Palmonari 9,10, Beatrice Panico 6, Alexandra Parmentier 1, Francesco Perfetto 6, Piergiorgio Picozza 1,2, Mirko Piersanti 11, Michele Pozzato 10, Ester Ricci 3,4, Marco Ricci 12, Sergio Bruno Ricciarini 7, Zouleikha Sahnoun 10, Valentina Scotti 6,13, Alessandro Sotgiu 1, Vincenzo Vitale 1, Simona Zoffoli <sup>14</sup> and Paolo Zuccon 3,4**


**Abstract:** The radiation belts in the Earth's magnetosphere pose a hazard to satellite systems and spacecraft missions (both manned and unmanned), heavily affecting payload design and resources, thus resulting in an impact on the overall mission performance and final costs. The NASA AE9/AP9/SPM radiation models for energetic electrons, protons, and plasma provide useful information on the near-Earth environment, but they are still incomplete as to some features and, for some energy ranges, their predictions are not based on a statistically sufficient sample of direct measurements. Therefore, it is of the upmost importance to provide new data and direct measurements to improve their output. In this work, the AP9 model is applied to the China Seismo-Electromagnetic Satellite (CSES-01) orbit to estimate the flux of energetic protons over the South Atlantic Anomaly during a short testing period of one day, 1 January 2021. Moreover, a preliminary comparison with

Bartocci, S.; Battiston, R.; Burger, W.J.; Campana, D.; Carfora, L.; Castellini, G.; Conti, L.; Contin, A.; et al. Trapped Proton Fluxes Estimation Inside the South Atlantic Anomaly Using the NASA AE9/AP9/SPM Radiation Models along the China Seismo-Electromagnetic Satellite Orbit. *Appl. Sci.* **2021**, *11*, 3465. https://doi.org/10.3390/app11083465

**Citation:** Martucci, M.; Sparvoli, R.;

Received: 9 March 2021 Accepted: 8 April 2021 Published: 13 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

proton data obtained from the High-Energy Particle Detector (HEPD) on board CSES-01 is carried out. This estimation will serve as the starting ground for a forthcoming complete data analysis, enabling extensive testing and validation of current theoretical and empirical models.

**Keywords:** trapped particles; South Atlantic Anomaly; AE9/AP9/SPM models; radiation belts

#### **1. Introduction**

The radiation belts, also known as Van Allen belts, are regions of the Earth's magnetosphere where energetic charged particles are subject to long-term magnetic trapping. The outer belt is mostly populated by electrons with hundreds of keV to MeV energies, while the inner belt mostly consists of an intense radiation of energetic protons (from MeV up to a few GeV), electrons/positrons (up to ∼8 MeV), and a minor component of ions [1,2]. Proton populations with energies above a few tens of MeV originate from the *β*-decay of free neutrons produced in the interaction between galactic cosmic-rays and the Earth's atmosphere in a mechanism called *Cosmic Ray Albedo Neutron Decay* (CRAND) [3,4]. Since the discovery of the Van Allen radiation belts, after the launch of the first Explorer satellites in 1958 and the Pioneer in 1959 [5,6], the scientific community has been considerably involved in modeling this space radiation environment. All these efforts were mostly aimed to meet the practical need of better understanding the significant radiation hazard to spacecraft and human crews. Several studies reported a direct association between the dynamic radiation environment and system or subsystem performances [7,8]. To address and solve these problems, more accurate, comprehensive, and up-to-date space radiation environment models have been developed by the National Reconnaissance Office (NRO) and the Air Force Research Laboratory (AFRL); the new AE9/AP9/SPM set of models for high-energy electrons, protons, and space plasma, respectively, is derived from 45 datasets obtained from sensors on board various satellites. These datasets have been processed to create maps of the particle fluxes along with estimates of uncertainties from both imperfect measurements and space weather variability [9]. A detailed comparison between the older AE8/AP8 and the newer AE9/AP9 models is reported in [10].

Gradual deterioration of spacecraft systems and components—and their overall performances—with accumulated dose is a fact, and various failures, due to phenomena associated with Single Event Effects (SEEs) or electrostatic discharge, are particularly common.

The first empirical models of the radiation belts, developed in the 1960s and 1970s by NASA, tried to describe and represent the radiation environment and their early versions, namely, AE8 and AP8 [11], were widely employed in spite of their limitations, especially at low altitudes [11–13]. Despite being successful in describing the radiation environment, even AE9/AP9/SPM are partly incomplete, and often their predictions are not based on a statistically sufficient sample of direct measurements [10]. For this reason, it is of key importance to test them and, above all, to provide new and reliable datasets from in-flight instruments to improve their output and accuracy. Among the scientific payloads on board the China Seismo-Electromagnetic Satellite (CSES-01), in Low-Earth Orbit since February 2018, the High-Energy Particle Detector (HEPD) has gone through an intense period of testing and calibration, and it is able to measure >3 MeV electrons and >35 MeV protons with high efficiency. With an overall expected mission duration of >5 years, and together with other similar missions planned in the coming years, measurements from HEPD could enable the testing and validation of the aforementioned models. In this work, we have used orbital information from CSES-01 and the AP9 model to estimate the flux of trapped protons over the South Atlantic Anomaly (SAA) during a short testing period, i.e., January 1, 2021, in order to assess the radiation level at CSES orbit in view of a comparison to experimental HEPD data in a forthcoming publication. A brief description of the South Atlantic Anomaly is given in Section 2, while some details on the CSES mission and the

HEPD payload are given in Section 3. The analysis is described in Section 4, results are presented in Section 5, and, finally, a brief discussion is presented in Section 6.

#### **2. The South Atlantic Anomaly**

The South Atlantic Anomaly (SAA) is one of the most well-known features of the Earth's magnetic field. It emerges as a consequence of the tilt (∼10◦) between the magnetic dipole axis of the Earth and its rotational axis and of the offset (∼500 km) between the dipole and the Earth centers. It can be considered as the response of an inverse flux path at the core–mantle boundary of the radial component of the geomagnetic field located approximately under the South Atlantic Ocean, which generates the hemisphere asymmetry of the Earth's magnetic field [14]; this region is characterized by an extremely low intensity of the geomagnetic field, and its behavior suggests that this asymmetry could be connected to the general decrease of the dipolar field and to the significant increase of the non-dipolar field in the Southern Atlantic region [15,16]. The extent area of the SAA at the surface of the Earth has been continuously growing since instrumental intensity measurements were made available. Several studies relate this as an indicator of a possible upcoming geomagnetic transition (excursion or reversal). It is generally accepted that such transitions are anticipated by flux patches of reversed polarity, slowly appearing at low- or mid-latitude, which migrate towards the pole [17,18].

The spatial and temporal evolution of the geomagnetic field has been monitored since 1832, when Carl-Friedrich Gauss performed the first intensity measurements in this region [19]. It has been shown that the magnetic dipole strength has been continuously decreasing [20], and data from the Swarm mission [21] revealed that two different patches are present over South America and near the coast of Africa, the latter growing at a rate of −2.54 × <sup>10</sup><sup>5</sup> nT per century [22]. A correct modelization of the SAA is of capital importance due to the high impact it has on human health and on instrumental efficiency [23]. Furthermore, recent studies indicate that the extent of the anomaly follows a log-periodic acceleration, resembling the behavior of a system that moves toward a critical transition [24].

#### **3. The CSES Scientific Mission**

The China Seismo-Electromagnetic Satellite (CSES-01) [25] is the first of a series of multi-instrument monitoring satellites scheduled for launch in the next few years; it is designed to study the near-Earth environment, addressing variations of the electromagnetic field, plasma parameters, and particle fluxes linked to natural sources or artificial emitters. The main scientific objective of this mission—resulting from a Chinese/Italian joint effort is to investigate possible correlations between the aforementioned perturbations and the occurrence of high-magnitude seismic events, but it is also well suited for studying a wide variety of space-weather phenomena triggered by solar–terrestrial interactions on short (i.e., geomagnetic storms, solar particle events, etc.) and long time-scales (i.e., cosmic-ray propagation, composition, etc.) [26,27]. A recent perspective [28] explained that claims based on self-organized criticality stating that at any moment any small earthquake can eventually cascade to a large event do not stand in view of the results obtained by natural time analysis [29,30]. The CSES-01 satellite, based on the Chinese 3-axis-stabilized CAST2000 platform (total mass ∼ 700 kg), is flying on a sun-synchronous polar orbit at a ∼507 km altitude, with a 97◦ inclination, a period of 94.6 min, and a 5-day revisiting periodicity. Nine scientific payloads are present on board CSES: two sets of particle detectors, namely, the High-Energy Particle Package (HEPP) [31] and the High-Energy Particle Detector (HEPD) [32]; a High-Precision Magnetometer (HPM) [33]; a Search-Coil Magnetometer (SCM) [34]; an Electric Field Detector (EFD) [35]; a Global Navigation Satellite System (GNSS) Occultation Receiver [36]; a Langmuir Probe (LAP) [37]; a Tre-Band Beacon transmitter (TBB) [38]; and a Plasma Analyzer Package (PAP) [39].

#### *The High-Energy Particle Detector*

The High-Energy Particle Detector (HEPD) is a light and compact (40.36 cm × 53.00 cm × 38.15 cm, total mass ∼45 kg) payload designed and built by the Limadou team, the Italian branch of the CSES Collaboration. The apparatus is made up of a series of sub-detectors:


The instrument is optimized to detect electrons in the 3 to 100 MeV energy range and protons between 35 and 250 MeV, as well as light nuclei. In these three years of flight, after a long period of calibration and testing, HEPD has been able to measure fluxes of low-energy galactic protons with great precision [40] and to observe the effects of the geomagnetic storm of August 2018 [41]. All the capabilities assessed in these years make HEPD well suited for the analysis of low-energy electrons and protons with good angular resolution and stability over time, which is particularly useful in highly anisotropic flux conditions like the ones encountered in SAA. More technical details can be found in [32,42,43].

#### **4. Materials and Methods**

The AE9/AP9/SPM set of models (version V1.50.001-release date December 2017) was downloaded from the Virtual Distributed Laboratory (VDL) website of the Air Force Research Laboratory (https://www.vdl.afrl.af.mil/programs/ae9ap9/, accessed on 1 February 2021). Element Set (ELSET) data—including Two-Line Elements (TLE) for the CSES satellite on 1 January 2021—have been retrieved from the Space-Track website (https://www.space-track.org/, accessed on 1 February 2021) and inserted in the code to generate the ephemeris of the satellite (ata5sresolution). The Simplified General Perturbations (*SGP4*) (the SGP4 propagator considers secular and periodic variations due the oblateness of the Earth, gravitational effects from Sun and Moon, and orbital drag and decay.) propagator has been preferred to the default *Kepler* with *J2* perturbation effects, for cross-checking purposes. Indeed, these orbital results have been further compared to the ones obtained using 2-min broadcast information downloaded from the satellite itself, to verify the correctness of the procedure. This cross-check includes the following:


After these steps, two sets of geographical/geomagnetic coordinates have been obtained: one calculated by NASA routines and another extracted by Limadou external routines used for trajectory propagation and magnetic field reconstruction since launch. A comparison between such datasets has been performed to assure the best possible agreement over the chosen testing period—January 1, 2021. The resulting discrepancies are <0.08% for LATs/LONs, <0.13% for altitude, and <0.21% for magnetic field intensity. Small discrepancies are probably due to the models using the last TLE entry for orbit

propagation, unlike external code picking the TLE closest in time. After ephemeris generation, an omnidirectional differential spectrum of trapped protons as a function of kinetic energy was created using AP9. This spectrum, averaged over all orbits of a single day, is shown in Figure 1; the blue arrow represents the HEPD low-energy thresholds for protons (30 MeV). Inside the inner radiation belt, trapped electron populations present a sharp threshold at ∼8 MeV, thus electrons of higher energy are virtually nonexistent. This means that in our future analysis, trapped protons will not be affected by any low-energy electrons contamination inside the SAA, consequently improving HEPD sensitivity to protons measurements.

**Figure 1.** Omnidirectional differential spectrum of trapped protons (blue squares) as a function of kinetic energy, obtained from the AP9 model and averaged over all the orbits of the testing day—1 January 2021. HEPD low-energy threshold for protons (30 MeV) is also depicted as a blue arrow.

The AP9 model could also be very useful to help define a fiducial area on the Earth's surface (longitude vs. latitude) that may be applied to the future HEPD data analysis of trapped particles. However, the geographical extension of the inner belt (and consequently of the SAA) is largely dependent on energy, as shown in Figure 2. The surface contours for the >1, >10, and >100 MeV trapped protons are depicted as red, blue, and green curves, respectively, in the panel. These contours highlight how the low-energy trapped proton component is distributed in the southern regions of the SAA-even superimposing to the outer belt, while higher energy populations are more clearly enclosed in the classical boundaries of the inner belt, i.e., in the area above Brazil and the Atlantic Ocean.

**Figure 2.** Geographical extension of the SAA for >1, >10, and >100 MeV protons (respectively, blue, red, and green curves in the panel), obtained from the AP9 model.

For a further, more precise comparison with experimental HEPD data, four different 20 min orbit portions were selected among those crossing the SAA , in order to build the related time-profiles of protons at various energies (With a period of 94.6 min, the satellite makes ~15 complete orbits per day.):


These passages over the SAA are represented in Figure 3 as a function of the geographical coordinates; while orbit 1 crosses the SAA in the outermost and peripheral region, orbit 3 crosses the bulk of the Anomaly, where particle fluxes are expected to be higher.

**Figure 3.** Representation of the four orbits chosen for the time-profile evaluation as a function of geographical latitude and longitude; orbit 1 appears to be more peripheral with respect to, for example, orbit 3.

#### **5. Results**

The differential, omnidirectional energy spectra of trapped protons along the four portions of orbits depicted in Figure 3, all generated by the AP9 model, are shown in Figure 4. As energy increases, the spectrum in each orbit decreases with a somewhat

different steepness, as expected. This is, in fact, due to the different aspects of the trapping mechanism, which is the resulting effects of cosmic ray albedo neutron decay (CRAND), solar proton injection, and radial diffusion [45]. The CRAND mechanism is the principal trapped proton source above ∼100 MeV, and the shape of the albedo neutron vertical spectrum above the geomagnetic cut-off is very similar to the one observed in trapped protons, i.e., the spectrum is decreasing as energy increases. On the contrary, the solar injection is more relevant below ~100 MeV (and it is more important for L > 2), while the radial diffusion tends to redistribute trapped particles in different L, so its effects are more complex. As a result, the trapped proton flux is strongly anisotropic and the overall spectrum changes rapidly, heavily depending on the region (latitude, longitude, L, etc.) where it is estimated.

**Figure 4.** Differential, omnidirectional energy spectra of >10 MeV trapped protons obtained with the AP9 model and averaged over each CSES-01 orbit (see title above each panel). In each orbit, the spectrum decreases as energy increases, as expected.

The time profiles (5-s resolution) for the >10 MeV trapped protons along the orbits defined above are also shown in Figure 5. In each panel, the color palette relates to the different particle energy. Note that during each portion of the CSES-01 orbits, the spectra possess a different shape as a function of time. For example, during orbit 1, trapped fluxes tend to decrease very rapidly, while for the other orbits this decrease is slower. This is due to the fact that orbit 1 crosses the SAA region only in its peripheral section, so the trapped population is only encountered for a small amount of time. Furthermore, in each panel it seems that energetic protons are more concentrated in the internal sectors of the SAA, while low-energy protons are more widely distributed and spread over a larger area; this was inferred also from Figure 2, and it is another proof of the high variability of trapped fluxes inside the inner belt.

**Figure 5.** Time profiles (5-s resolution) of 10–300 MeV trapped protons estimated from the AP9 model along CSES orbits 1–4 (from the top left panel to the lower right panel). As expected, higher energies have lower fluxes, while lower energies tend to have higher fluxes.

To assess the level of agreement between the AP9 model and the experimental data, a preliminary analysis was conducted using omnidirectional ∼50 MeV calibrated proton data from HEPD, obtained following the same procedure used in [40]. As can be seen from the four panels in Figure 6, the agreement seems to be good, even if some small discrepancies are evident, mostly in the peripheral regions of the SAA; these are probably due to the different operational definition of South Atlantic Anomaly that was used to derive the data with HEPD (For the HEPD data analysis, we define the South Atlantic Anomaly as the region enclosed in a value of the magnetic field >20,000 nT.). Further studies are needed to verify the agreement even in a longer time period and with the extensive use of simulations.

For the purposes of this preliminary analysis, the uncertainties of the AP9 model related to measurement, gap-filling, or dynamic variations due to space-weather processes are not taken into account. A future, more complete comparison with HEPD observations will require a precise assessment of the AP9 confidence levels, in order to better evaluate the match with experimental data. Considering that CSES-01 will be operative in a period of strong minimum between the end of the 24th solar cycle and the start of the 25th, no major effect related to space weather variability is expected.

**Figure 6.** Time profiles (5-s resolution) of 50 MeV trapped protons estimated from the Ap9 model and compared with preliminary data of ∼50 MeV proton data (black circles) from the HEPD instrument on board the CSES-01 satellite. The analysis has been carried out using the procedure described in [40]. The agreement between the data and the model appears generally good, despite showing small discrepancies, especially in the peripheral regions of the SAA. Only statistical errors are reported.

#### **6. Discussion and Conclusions**

The NASA AE9/AP9/SPM set of models represents an important approach to specify the radiation environment for modern satellite design applications. In this work, this suite of models has been employed to estimate trapped proton fluxes over the South Atlantic Anomaly for some orbits of the China Seismo-Electromagnetic Satellite on 1 January 2021. This is intended as a starting point for a future analysis that will include data from the High-Energy Particle Detector. After three years of calibration and testing, HEPD has proven capable of measuring low-energy particle fluxes (>3 MeV electrons and >35 MeV protons) with precision and stability over time; among the others, these two characteristics in particular are very suitable for the measurement of strongly anisotropic particle fluxes, such as those trapped in SAA. Thus, HEPD, together with the other payloads on board CSES (such as those of the HEPP suite), can provide excellent cross-calibration for these radiation environment models at LEO. A preliminary analysis on HEPD proton data has been conducted to assess the agreement between the AP9 model and experimental data, and it seems already acceptable, even if some discrepancies—that need to be studied—are present. It is important to remember, as already mentioned, that there is also a certain number of known issues in these models:


Moreover, much of the validation of these models was performed using the Van Allen Probe mission [46], which provided a rich set of energetic particle and plasma data from the many instruments the spacecraft carried on board, together with a good pitch angle and energy resolution; unfortunately, after the end of the mission, new data are necessary to continue validation and to explore higher energy ranges with more statistics. HEPD proved to be able to cover this role, performing measurements with precision and stability in time; besides, new CSES missions (with more HEPD-like particle detectors) are already planned for the next years, greatly expanding the data-taking period by several years into the 2020s.

**Author Contributions:** Writing original draft, M.M. (Matteo Martucci); Conceptualization, R.S; Writing-review and editing, R.S., R.B., L.C. (Livio Conti), F.M.F., M.P. (Mirko Piersanti), A.P., C.D.D., A.S., F.P. (Francesco Palma); Designing the experiment or calibration or data production and processing, S.B., R.B., W.J.B., D.C., L.C. (Luca Carfora), G.C., L.C. (Livio Conti), A.C., C.D.D., C.D.S., F.M.F., R.I., I.L., N.M., G.M., M.M. (Matteo Mergé), A.O., G.O., F.P. (Francesco Palma), F.P. (Federico Palmonari), B.P., A.P., F.P. (Francesco Perfetto), P.P., M.P. (Mirko Piersanti), M.P. (Michele Pozzato), E.R., M.R., S.B.R., Z.S., V.S., A.S., V.V., S.Z., P.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** CSES/HEPD data can be found in www.leos.ac.cn/, accessed on 1 February 2021.

**Acknowledgments:** This work makes use of data from the CSES mission (www.leos.ac.cn/, accessed on 1 February 2021), a project funded by China National Space Administration (CNSA), China Earthquake Administration (CEA) in collaboration with the Italian Space Agency (ASI), National Institute for Nuclear Physics (INFN), Institute for Applied Physics (IFAC-CNR), and Institute for Space Astrophysics and Planetology (INAF-IAPS). We kindly acknowledge the AFRL for providing the AE9/Ap9/SPM set of models. This work was supported by the Italian Space Agency in the framework of the "Accordo Attuativo 2020-32.HH.0 Limadou Scienza+" (CUP F19C20000110005) and the ASI-INFN agreement n.2014-037-R.0, addendum 2014-037-R-1-2017.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


### *Article* **Low-Pass Filtering Method for Poisson Data Time Series**

**Victor Getmanov 1,2,\*, Vladislav Chinkin 1, Roman Sidorov 1,\*, Alexei Gvishiani 1,2, Mikhail Dobrovolsky 1, Anatoly Soloviev 1,2, Anna Dmitrieva 1,3, Anna Kovylyaeva 1,3, Nataliya Osetrova 1,3 and Igor Yashin 1,3**


**Abstract:** Problems of digital processing of Poisson-distributed data time series from various counters of radiation particles, photons, slow neutrons etc. are relevant for experimental physics and measuring technology. A low-pass filtering method for normalized Poisson-distributed data time series is proposed. A digital quasi-Gaussian filter is designed, with a finite impulse response and non-negative weights. The quasi-Gaussian filter synthesis is implemented using the technology of stochastic global minimization and modification of the annealing simulation algorithm. The results of testing the filtering method and the quasi-Gaussian filter on model and experimental normalized Poisson data from the URAGAN muon hodoscope, that have confirmed their effectiveness, are presented.

**Keywords:** Poisson data; time series; quasi-Gaussian filter; digital filtering; optimization; global minimization; annealing simulation algorithm

#### **1. Introduction**

The article proposes a low-pass filtering method for Poisson-distributed data time series, based on a specially developed digital low-pass filter with finite impulse response (FIR filter), with gain equal to one at zero frequencies and non-negative weighting factors.

Here, low-pass filtering is applied in order to reduce noise in Poisson-distributed data to ensure the recognition of emerging fluctuations of mathematical expectations in them. Poisson-distributed, or Poisson data are found in various physical systems, for example, related to the heliosphere and magnetosphere of the Earth; the fluctuations of mathematical expectations of these data may contain information regarding the structures and characteristics of these systems.

A particular feature of the Poisson data origin is that they contain sufficient noises; it is known, for example, from [1] that their variance is numerically equal to mathematical expectation. Noise reduction in Poisson data can be achieved using common FIR filters [2,3], to which, within the framework of this article, we refer the filters based on commonly used windowing techniques, frequency sampling and inverse Fourier transforms [4,5]. However, there are a number of scientific and technical problems for which their application is not fully effective, for example, (1) recognition of small (in size and duration) mathematical expectation fluctuations in Poisson datasets; (2) digital processing of Poisson data with small mathematical expectation values.

Common FIR filters can potentially be used for the mentioned tasks, and their synthesis can be implemented according to given dimensions and cutoff frequencies. The synthesis procedures for common FIR filters are, in essence, the variants of approximation

**Citation:** Getmanov, V.; Chinkin, V.; Sidorov, R.; Gvishiani, A.; Dobrovolsky, M.; Soloviev, A.; Dmitrieva, A.; Kovylyaeva, A.; Osetrova, N.; Yashin, I. Low-Pass Filtering Method for Poisson Data Time Series. *Appl. Sci.* **2021**, *11*, 4524. https://doi.org/10.3390/app11104524

Academic Editors: Roberta Sparvoli and Matteo Martucci

Received: 23 April 2021 Accepted: 12 May 2021 Published: 15 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

procedures for the specified species frequency response (FR) types; the accuracy of the FR approximations depends on the specified dimensions for the synthesized filters. Obviously, at large dimensions, the accuracy of these approximations is high and the errors in the resulting cutoff frequencies are small. For the case of small dimensions, the approximation accuracy turns out to be low and, as a consequence, cutoff frequencies are realized with significant errors which prevent low-pass filtering. We can assume that the filtering procedure proposed here should be performed by filters with low dimensions and cutoff frequencies and with gain values equal to one in order to avoid mathematical expectation distortions, and with non-negative weight factors in order to provide non-negativity of filtering results taking into account the Poisson property of the data.

The indicated problem leads to the need to formulate the synthesis problem for a special digital low-pass FIR filter, which takes into account the requirements—restrictions on dimensionality, cutoff frequency, gain at zero frequencies, and weighting factors.

Here, a FIR filter is proposed, which is further denoted as a quasi-Gaussian filter, the frequency response of which is formed on the basis of approximating a Gaussian function and ensuring the implementation of the mentioned constraints conditions using a special optimization method.

Gaussian filters, the frequency response of which is implemented based on the approximation of the Gaussian function, are widely used in modern scientific and technical practice [6,7]. However, as a rule, the known variants of Gaussian filters with the approximation of the frequency response do not take into account the above-mentioned conditions (restrictions).

Problems of digital processing of Poisson data time series from muon counters in muon detectors and telescopes [8], counters of elementary particles of alpha-beta-gamma radiation, photon counters, slow neutrons, etc. [9], taking into account their specificity, are relevant for experimental physics. Digital processing of Poisson data, including the Gaussian filtering application, can be outside of experimental physics, for example, in medical technology for imaging blood vessels and tumor therapy with particle beams, in measuring technology for tribological studies of the surfaces of metal parts, in astronomy for gamma telescopes, in muon tomography for recognizing cavities in rocks, and building structures and many other applications.

One of the applications of the designed filter proposed here is the digital processing of the data from the URAGAN muon hodoscope (MH) designed by NRNU MEPhI [10,11]. The MH is a computerized measuring device that estimates the intensities of muon fluxes by counting the number of elementary particles—muons—registered in its detector for a set of solid angles with a set time step. Within the framework of this article, MH can be interpreted as a distributed set of muon counters, consisting of primary and secondary information converters.

From each primary MH transducer, the initial Poisson data—time series of random non-negative integers *N*(*Tk*)—the quantities of Poisson-distributed events recorded in a given solid angle at time intervals (*T*(*k* − 1), *T*(*k* − 1) + *T*0*k*), *k* = 1, 2, ... , *k*0, where *T* = 1 minute. Due to the features of the MH design, registration intervals *T*0*k*are random with a uniform distribution law in the range *T*0 min ≤ *T*0*<sup>k</sup>* ≤ *T*0 max < *T*.

From each secondary MH transducer, the 1-minute-sampled normalized Poisson data *Y*(*Tk*) are generated for a given solid angle by reducing to one second and calculating the averaged normalized Poisson data *Y*(*T*0*n*) with an hourly discreteness according to the following relations:

$$Y(Tk) = N(Tk) / T\_{0k}, \\ Y(T\_0 n) = \frac{1}{60} \sum\_{k=1+60(n-1)}^{k=60n} Y(Tk), \\ n = 1, 2, \dots, T\_0 = T \cdot 60. \tag{1}$$

Data resulting from (1) are produced for the whole set of solid angles; next, they are placed into time series of matrix MH data in the database [12].

#### **2. Method**

*2.1. Quasi–Gaussian Digital Low-Pass Filter*

2.1.1. Statement of the Problem

One-dimensional FIR filter synthesized here is built according to the following difference equation:

$$X(T\_0 n) = \sum\_{s=0}^{s\_0} a\_s Y(T\_0(n-s)), n = 1, 2, \dots,\tag{2}$$

where *r*<sup>0</sup> = *s*<sup>0</sup> + 1 is the FIR filter dimension, *a<sup>T</sup>* = (*a*0, *a*1, ... , *as*<sup>0</sup> ) is a weight factors vector, *X*(*T*0*n*) is the output time series, *Y*(*T*0*n*) is the FIR filter input—the hourly normalized Poisson data time series from MH according to (2), which begins from the values *Y*(*T*0(1 − *s*0)), *Y*(*T*0(1− *s*<sup>0</sup> +1)), *Y*(*T*0(1− *s*<sup>0</sup> +2)), .... Transfer function (TF) *H*(*jωT*0, *a*) for filter(2) is defined as follows:

$$H(j\omega T\_{0\prime}a) = \sum\_{s=0}^{s\_0} a\_s e^{-j2\pi\omega \cdot T\_{0\prime}}.\tag{3}$$

Here *ω* is the TF frequency parameter. For (3), a normalized fequence is introduced, *w*, *ωT*<sup>0</sup> = *wπ*, 0 ≤ *w* ≤ 1.0, and its discrete values are calculated: *wl*

$$dw = 1.0/L\_0, \ w\_l = dw(l-1), \ l = 1, \dots, L\_l \\ L = L\_0 + 1. \tag{4}$$

The frequency response (FR) *H*(*wl*, *a*) = |*H*(*jwl*, *a*)|, considering (3), is the following:

$$\begin{aligned} H(w\_{l\prime}, a)^2 &= H\_1^2(w\_{l\prime}a) + H\_2^2(w\_{l\prime}a), \\ H\_1(w\_{l\prime}a) &= \sum\_{s=0}^{s\_0} a\_s \cos(2\pi \tau w\_{l\prime}s), \\ H\_2(w\_{l\prime}a) &= \sum\_{s=0}^{s\_0} a\_s \sin(2\pi \tau w\_{l\prime}s) \end{aligned} \tag{5}$$

for discrete normalized frequencies *wl*, *l* = 1, ... , *L* according to (4). The cutoff frequency *wc* for FR (5) is found based on the equality |*H*(*jwA*, *a*)| <sup>2</sup> = 0.5.

For a low-frequency FIR filter synthesis, the FR of the prototype filter is used, based on a Gaussian function *H*0*g*(*w*, *wc*0)

$$H\_{0\lg}(w, w\_{c0}) = \exp(-(w/w\_{c0})^2). \tag{6}$$

2.1.2. Synthesis Requiements

The problem of synthesis of the supposed FIR filter is solved based on the approximation of the FR function *H*0*g*(*wl*, *wc*) (6) in discrete points *wl*, *l* = 1, ... , *L* with a FR function *Hg*(*wl*, 0) according to (5). A functional *S*(*H*0*g*, *a*, *wc*) is formed:

$$S(H\_{\mathbb{Q}\mathcal{Y}}, a, w\_{\varepsilon}) = \sum\_{l=1}^{L} [(\sum\_{s=0}^{\mathbb{S}\_{0}} a\_{s} \mathbb{C}\_{s}(w\_{l}))^{2} + (\sum\_{s=0}^{\mathbb{S}\_{0}} a\_{s} \mathbb{S}\_{s}(w\_{l}))^{2} - H\_{\mathbb{Q}\mathcal{Y}}^{2}(w\_{l}, w\_{\varepsilon})]^{2}.\tag{7}$$

Obviously, the FR (5) represents a function which is polyharmonic in frequency *wl*. In case the prototype filter FR frequency derivative has discontinuities or is subject to strong alternations, e.g., if FR is a trapezoidal function, then the FR of the synthesized FIR filter, obtained based on approximation, will contain fluctuations due to the so-called Gibbs effect. Elimination and reduction of these fluctuations are usually achieved by choosing a suitable smooth prototype filter FR function. The smoothness requirement is largely satisfied by the Gaussian function (6).It should be noted that the Gaussian function is naturally suitable for the FR of a low-pass filter, since its values (6) practically differ from zero only in the region of low frequencies.

The requirements listed in the Introduction lead to formalized requirements:

a. Ensuring that the gain at zero frequencies is equal to one:

$$A = H(0, a) = \sum\_{s=0}^{80} a\_s, a \in A\_1, \ A\_1 = \{a : (1 = \sum\_{s=0}^{80} a\_s)\};\tag{8}$$

b. Ensuring non-negativity of coefficients:

$$a \in A\_{0\prime} \quad A\_0 = \{ a : (0 \le a\_s, s = 0, 1, \dots, s\_0) \};\tag{9}$$

For the synthesis procedure, it is assumed to set a small value *r*0, based on the a priori known duration of fluctuations, and some small cutoff frequency value *wc* for a prototype filter. The quasi-Gaussian filter synthesis procedure, consisting of finding the optimal coefficients *a*◦ *<sup>s</sup>* , *s* = 0, 1, ... ,*s*0, taking into account the requirements **a**,**b**, Equations (8) and (9) the predefined *r*0and *wc*, is performed on the basis of the approximation problem, which reduces to the implementation of conditional minimization:

$$a^\odot(w\_\mathfrak{c}) = \arg \{ \min\_{a \in A\_0, a \in A\_1} S(H\_{0 \lg \prime} a, w\_\mathfrak{c}) \}. \tag{10}$$

For a given small dimension *r*<sup>0</sup> of the synthesized quasi-Gaussian filter and a given small cutoff frequency *wc* for a prototype filter, the value for cutoff frequency to be found for a quasi-Gaussian filter is *wcg*, and the filter FR for the frequencies *wl* is denoted as *Hg*(*wl*, *wcg*, *a*◦), *l* = 1, . . . , *L* .

The minimization of (10) could be performed based on modified direct zero-order optimization methods, taking into account the restrictions (8) and (9). However, because the (7) functional is multi-extremal, traditional modified direct methods, for example, using the coordinate descent method, the Hook–Jeeves method, the random descent method, etc. [13] do not provide successful minimization. The listed methods, as a rule, lead to "getting stuck" with search procedures in local minima.

#### *2.2. Quasi–Gaussian Filter Synthesis Procedure*

We can synthesize the quasi-Gaussian filter based on the technology of stochastic global minimization of the (7) functional with the constraints (8) and (9) using the optimization algorithm for annealing simulation [14,15]. To implement it, we will use the simulannealbnd.mat software module from the Matlab Global Optimization Toolbox [16].

Let us form a parallelepiped of constraints *A*<sup>0</sup> of dimension *r*<sup>0</sup> with boundaries *ar*, *r* = 1, ... ,*r*0—*a* ∈ *A*0, *A*<sup>0</sup> = {*a* : (0 ≤ *ar* ≤ *ar*, *r* = 1, ... ,*r*0)} and a new—with respect to (7)—functional *S*(*H*0*g*, *a*) with a penalty term taking into account the constraint equality (8). Let us implement the global minimization of *S*(*H*0*g*, *a*) taking into account *A*<sup>0</sup> using [16].

Let us set the initial vectors for the first iteration *a*1(*I*) ∈ *A*0, uniformly distributed in *A*0, *I*—a single descent procedure , *I* = 1, 2, ... , *I*0, *I*0—a total number of descent procedures. Let us assume that each descent procedure consists of *m*0—a total number of iterations, *m*—a single iteration, *m* = 1, 2, ... , *m*0. During descent, we assume that the initial value of the vector of parameters for (*m* + 1)-st iteration is equal to the calculated optimal value for the vector of parameters for *m*-th iteration—*am*+1(*I*) = *a*◦ *<sup>m</sup>*(*I*) . In each iteration, we perform *n*<sup>0</sup> descent steps, *n* is a descent step number, *n* = 1, 2, ... , *n*0. Next, we will calculate the sequence of the functional *S*(*m*0, *I*)= *S*(*H*0*g*, *a*◦) values and the corresponding optimal vectors *a*◦(*m*0, *I*), *I* = 1, 2, ... , *I*0. For global minimization, we search for the optimal index *I*◦ corresponding to the minimum of the *S*(*m*0, *I*) functional, and the optimal vector *a*◦ using brute force:

$$I^\diamond = \arg \left( \min\_{1 \le I \le I\_0} S(m\_{0\prime}I) \right), \\ a^\diamond = a^\diamond (m\_{0\prime}I^\diamond).$$

On Figure 1, the example plots of the minimized *S*(*m*, *I*) functionals are displayed, depending on iteration number *m* and the descent procedure number *I*. Functionals are

shown starting with *m* = 2, since for *m* = 1 their values are very large. Here, *m*<sup>0</sup> = 20; as the iteration number increases, the values of the functionals decrease. During the optimization process, a movement is made in a *r*0-dimensional space from one local minimum to another.

**Figure 1.** Plots of descent procedures—minimization of functionals *S*(*m*, *I*), *I* = 1, 2, ... , *I*0, *m* = 1, 2, . . . , *m*0.

Let us consider an example of quasi-Gaussian FIR filter synthesis. Based on the analysis of hourly experimental MH data from [12], it was found that the durations of possible fluctuations of the mathematical expectation in them were, on average, ≈10 ÷ 20 h and more. The dimension value *r*0, that could possibly allow the recognition of such fluctuations in mathematical expectations, was equal to 8. For a prototype filter FR (6), the parameter *wc*<sup>0</sup> was related to the assigned cutoff frequency *wc* based on (6)

$$(0.5)^{1/2} = \exp(-(w\_{\varepsilon}/w\_{\varepsilon0})^2), \\ w\_{\varepsilon0} = w\_{\varepsilon}/(0.5 \cdot \ln 2)^{1/2}.$$

We assign the cutoff frequency *wc* = 0.1, find *wc*<sup>0</sup> and define *H*0*g*(*w*, *wc*)—the prototype filter FR. By defining *L* we set the number of discrete normalized frequencies *wl* of calculations of the functional (7) for 0 ≤ *wl* ≤ 1.0, let us assume that *L* = 100 in our calculations. The polyharmonic FR function |*H*(*jw*)| (5) is formed from components performing 1, 2, . . . , *s*<sup>0</sup> fluctuations in this interval. For the accepted values *L* and *r*0, one period of the polyharmonic component with the maximum frequency corresponding to the number *s*<sup>0</sup> in (5), accounted for ≈15 sampling points of normalized frequencies *wl*, *l* = 1, ... , *L*, which fully provided a fairly accurate calculation of the functional (7) necessary for direct search.

Let us calculate the vector of factors *a*◦, form the synthesized quasi-Gaussian filter FR *Hg*(*w*, *wcg*, *a*◦) and define the cutoff frequency *wcg* = 0.175 .

For the comparison, let us synthesize a common FIR filter using the fir1.mat module [3]. For the dimension *r*<sup>0</sup> = 8 and the assigned cutoff frequency *wc* = 0.1 we find out the final cutoff frequency *wc f* = 0.275; let us denote the FR as *Hf*(*w*, *wc f*). On Figure 2, the FR plots for *H*0*g*(*w*, *wc*), *Hg*(*w*, *wcg*, *a*◦), *Hf*(*w*, *wc f*) are displayed. It is seen that, in case of low *r*<sup>0</sup> , the quasi-Gaussian filter FR was characterized by a better approximation to the prototype filter FR than the one of the common FIR filter.

**Figure 2.** FR plots: *H*0*g*(*w*, *wc*) (line 1), *Hg*(*w*, *wcg*, *a*◦) (line 2), *Hf*(*w*, *wc f*) (line 3).

Note that the proposed FIR filter, with the same dimension as the common FIR filter, made it possible to provide a lower value of the cutoff frequency than the realized cutoff frequency for the common FIR filter. The calculated cutoff frequencies of resulting FRs for common FIR filters synthesized using frequency sampling method and Fourier transforms [3,4] insignificantly (by ≈5–7% ) differ from the cutoff frequency *wc f* = 0.275. This gives a reason to make a conclusion about the advantages of a quasi-Gaussian filter over standard FIR filters.

#### **3. Results**

#### *3.1. Testing the Method and the Quasi–Gaussian Filter on Model Normalized Poisson Data*

#### 3.1.1. Testing on Model Hourly Data Using Statistical Modeling

Testing of the proposed method and quasi-Gaussian filter was carried out on model hourly normalized data using statistical modeling [17]. For this purpose, on the basis of the Matlab module exprnd.mat [18], exponentially distributed model random numbers *τi*, *i* = 1, 2, ... were generated, with their mean value *τM*<sup>0</sup> , and the evenly distributed random registration time intervals *T*0*k*, *k* = 1, 2, ... within the range *T*<sup>0</sup> *min* ≤ *T*0*<sup>k</sup>* ≤ *T*<sup>0</sup> *max*. The number of Poisson model events *NM*(*Tk*) was counted on the registration time intervals *T*0*k*. Finding *NM*(*Tk*) was carried out by solving the conditional maximization problems:

$$N\_M(Tk) = \arg\{\max N\_M\}\_i > 0,\tag{11}$$

providing that *<sup>T</sup>*0*<sup>k</sup>* <sup>−</sup> <sup>∑</sup>*NM <sup>i</sup>*=<sup>1</sup> *τ*, where for *T*0*k*, *k* = 1, 2, ... *k*<sup>0</sup> the range bounds *T*<sup>0</sup> *min* = 57 s, *T*<sup>0</sup> *max* = 59.5 s were assigned (see the Introduction section). Initial model 1-minutesampled and normalized Poisson-distributed data were constructed according to (11) and the calculation of relations *NM*(*Tk*), similar to (1):

$$
\overline{N}\_M(Tk) = N\_M(Tk) / T\_{0k}, \\
k = 1, 2, \dots, k\_0. \tag{12}
$$

The modulation of the average number of Poisson events in order to model decreases (increases) in the mathematical expectation was carried out by specifying the mean value function *τM*0(*Tk*) on the intervals (*T*(*k* − 1), *Tk*) for *k* from (12). For this, the relative modulation function *μ*(*Tk*), *k* = 1, 2, ... , *k*<sup>0</sup> was formed and the initial temporal index of the modulation decrease *ka*, the duration of the decrease *dka* and the depth of the relative decrease *dμ* ;. The function *μ*(*Tk*) was represented by the relations *μ*(*Tk*) = 1 − *dμ* for *ka* ≤ *k* ≤ *ka* + *dka*, *μ*(*Tk*) = 1 for 1 ≤ *k* < *ka*, *ka* + *dka* + 1 ≤ *k* ≤ *k*0.

For the calculation example, the average number of Poisson model events per minute was set *NM*<sup>0</sup> = 25, normalized average *NM*<sup>0</sup> = *NM*0/*T*, modulated normalized mean *NM*0(*Tk*) = *NM*0*μ*(*Tk*) = *NM*0*μ*(*Tk*)/*T*, *k* = 1, 2, ... , *k*<sup>0</sup> and the parameter *τM*0(*Tk*) = 1/(*NM*0(*Tk*) − 1) was calculated.

Based on [18], random exponentially distributed numbers with *τM*0(*Tk*) and random evenly distributed values with *T*<sup>0</sup> *min* = 57*s*, *T*<sup>0</sup> *max* = 59, 5*s* were generated, with the use of which by (11), model Poisson data *NM*(*Tk*) and by (12)—normalized Poisson data *NM*(*Tk*) were calculated. Further, similarly to (1), a time series of averaged model hourly normalized Poisson data was formed:

$$\mathcal{Y}\_M(T\_0 n) = \frac{1}{60} \sum\_{k=1+60(n-1)}^{k=60n} \overline{\mathcal{N}}\_M(Tk), n = 1, 2, \dots, n\_0, n\_0 = \text{ent}(k\_0/60). \tag{13}$$

For modeling, we assumed *k*<sup>0</sup> = 6000, which corresponded to the model minute data produced during 4.166 days. For the modulation function, the values *ka* = 1920, *dka* = 1440 and *dμ* = 0.02 were taken. Model hourly averaged data *YM*(*T*0*n*) for (13) with *n*<sup>0</sup> = 100, *na*<sup>1</sup> = 32, *na*<sup>2</sup> = *n*<sup>1</sup> + *dna dna* = 24.

Figure 3 shows an example of statistical modeling results: the jagged light gray line with index 1 displays the *YE*(*T*0*n*) plot; the solid line with index 2 denotes the fragment of *XEG*(*T*0*n*) which is the result of filtering the model dataset using a quasi-Gaussian filter; for comparison, the dashed line with index 3 denotes the fragment *XEF*(*T*0*n*) which is the result of filtering the model dataset using the software module fir1.mat [3]. Model piecewise constant modulating function *YM*0(*T*0*n*) = *NM*(*T*0*n*), represented by a dotted line (index 4), *m*<sup>0</sup> + *dm* = *YM*0(*T*0*n*) = 0.4165 for 1 ≤ *n* < *na*1, *na*<sup>2</sup> ≤ *n* < *n*0, *<sup>m</sup>*<sup>0</sup> = *YM*0(*T*0*n*) = 0.4087 for *na*<sup>1</sup> ≤ *<sup>n</sup>* ≤ *na*2, where the value of *dm* = 0.833 × <sup>10</sup>−<sup>2</sup> corresponded to the predefined 2% decrease. The plots show that the result of the quasi-Gaussian filter application (line 2) is a better approximation to the model piecewise constant modulation (line 4) than the result of a common FIR filtering (line 3).

**Figure 3.** Fragments of model datasets *YM*(*T*0*n*) (line 1), filtering results *XMG*(*T*0*n*) (line 2), *XMF*(*T*0*n*) (line 3) and model modulating function *YM*0(*T*0*n*) (line 4).

The calculation of approximate estimates of filtering errors for the quasi-Gaussian filter and fir1-filter was performed by calculating the root-mean-square (RMS) errors according to the following formulas for datasets *YM*0(*T*0*n*), *YMG*(*T*0*n*), *YMF*(*T*0*n*):

$$
\sigma\_{MG}^2 = \frac{1}{n\_0} \sum\_{n=-1}^{n\_0} \left( \chi\_{M0}(T\_0 n) - \chi\_{MG}(T\_0 n) \right)^2,\\
\sigma\_{MF}^2 = \frac{1}{n\_0} \sum\_{n=-1}^{n\_0} \left( \chi\_{M0}(T\_0 n) - \chi\_{MF}(T\_0 n) \right)^2. \tag{14}
$$

Results of a large number of tests performed for (14) showed that the *σMG* error values for *XMG*(*T*0*i*) regarding *YMG*(*T*0*n*) are, on average, 15–30% less than the corresponding *σMF* error values for *XMF*(*T*0*i*). An overview of model *XMG*(*T*0*i*) and *XMF*(*T*0*i*) (Figure 3) made it possible to ensure that the minimum duration of the interval, within which recognition for the decrease *dμ* = 0.02 can be performed, is 12–24 h.

The proposed method and the quasi-Gaussian filter provided more noise reduction than a common FIR filter. Consideration of the results of statistical modeling made it possible to draw a conclusion about the efficiency of the quasi-Gaussian filtering method.

#### 3.1.2. Estimation of Mathematical Expectation and Its Root Mean Square Errors

Testing of the method and quasi-Gaussian filter for estimating the mathematical expectation and the RMS of its errors depending on *dna*—the duration of decreases and *dμ*—the relative decrease value were carried out using statistical tests [17]. Random datasets *YM*(*s*, *T*0*i*), *XMG*(*s*, *T*0*i*), *XMF*(*s*, *T*0*i*) , *s* = 1, 2, ... , *M*, where *s* is the number of the dataset, *M* is the total quantity of datasets. The estimates of mathematical expectation *m*◦ *<sup>g</sup>*(*dna*, *dμ*) and RMS values *σ*◦ *<sup>g</sup>* (*dna*, *dμ*) for *XMG*(*s*, *T*0*i*) for a set of values *dna* and *dμ*

$$m\_{\mathcal{S}}^{\diamond}(s, dn\_{\mathfrak{a}}, d\mu) = \frac{1}{n\_{\mathfrak{a}}} \sum\_{n=-n\_{\mathfrak{a}}}^{n\_{\mathfrak{a}}+d n\_{\mathfrak{a}}} X\_{\mathrm{MG}}(s, T\_{0}n), m\_{\mathcal{S}}^{\diamond}(dn\_{\mathfrak{a}}, d\mu) \ = \frac{1}{M} \sum\_{s=1}^{M} m\_{\mathcal{S}}^{\diamond}(s, dn\_{\mathfrak{a}}, d\mu),$$

$$\sigma\_{\mathcal{S}}^{\diamond}(s, dn\_{\mathfrak{a}}, d\mu) \ = \frac{1}{n\_{\mathfrak{a}} - 1} \sum\_{n=-n\_{\mathfrak{a}}}^{n\_{\mathfrak{a}} + d n\_{\mathfrak{a}}} (X\_{\mathrm{MG}}(s, T\_{0}n) - m\_{\mathcal{S}}^{\diamond}(s, dn\_{\mathfrak{a}}, d\mu))^{2}, \quad \text{(15)}$$

$$\sigma\_{\mathcal{S}}^{\diamond}(d n\_{\mathfrak{a}}, d\mu) \ = \frac{1}{M} \sum\_{s=1}^{M} \sigma\_{\mathcal{S}}^{\diamond}(s, dn\_{\mathfrak{a}}, d\mu).$$

The coefficients of relative errors *ε*◦ *gm*(*dna*, *dμ*), *ε*◦ *<sup>g</sup>σ*(*dna*, *dμ*) of the quasi-Gaussian filter as ratios of errors *m*◦ *<sup>g</sup>*(*dna*, *dμ*) − *m*<sup>0</sup> and RMS *σ*◦ *<sup>g</sup>* (*s*, *dna*, *dμ*) to the values of *dm* reductions are the following:

$$\varepsilon^{\mathbb{C}}\_{\mathbb{S}^{m}}(dn\_{4},d\mu) = (m^{\mathbb{C}}\_{\mathbb{S}}(dn\_{4},d\mu) - m\_{0})/dm,\\ \varepsilon^{\mathbb{C}}\_{\mathbb{S}^{\mathbb{D}}}(dn\_{4},d\mu) = (\upsilon^{\mathbb{C}}\_{\mathbb{S}}(dn\_{4},d\mu))/dm\dots \tag{16}$$

The coefficients *ε*◦ *gm*, *ε*◦ *<sup>g</sup>σ*, calculated for *dna*, *dμ*, characterized the recognition capabilities of quasi-Gaussian filtering model decreases. Similarly, using (15) and (16) *m*◦ *<sup>f</sup>*(*dna*, *dμ*) and *σ*◦ *<sup>f</sup>* (*dna*, *dμ*) for *XMF*(*s*, *T*0*n*) and the coefficients *ε*◦ *f m*(*dna*, *dμ*), *ε*◦ *<sup>f</sup> <sup>σ</sup>*(*dna*, *dμ*). On Figure 4, the results of statistical tests are displayed, where *M* = 500. The *ε*◦ *gm*(*dna*, *dμ*) coefficients plots are the solid lines with indices 1, 2, and the *ε*◦ *f m*(*dna*, *dμ*) plots are the dashed lines with indices 3, 4. The coefficients *ε*◦ *gm*, *ε*◦ *f m* are given depending on the duration with the values *dna* = 12, 24, 48, 72 h and relative decreases in *dμ*, taking the values of 0.01 (indices 1, 3) and 0.03 (indices 2, 4).

The effect of quasi-Gaussian filtering was determined based on the calculation of *δε*◦ *f g*,*m*—the rates of errors with respect to the mathematical expectations:

$$\delta \varepsilon\_{f\S,m}^{\diamond}(d\mathfrak{n}\_a, d\mu) \;= \left(\varepsilon\_{fm}^{\diamond}(d\mathfrak{n}\_a, d\mu) - \varepsilon\_{\mathfrak{g}m}^{\diamond}(d\mathfrak{n}\_a, d\mu)\right) / \varepsilon\_{\mathfrak{g}m}^{\diamond}(d\mathfrak{n}\_a, d\mu) \tag{17}$$

The results of the *δε*◦ *f g*,*<sup>m</sup>* calculations according to (17) for some *dμ* and *dna* values are:


Analysis of the error values showed that the *ε*◦ *gm* rate values appeared to be about 10–30% lower than the *ε*◦ *f m* values. The nature of the dependencies of the estimates of the error coefficients for the *ε*◦ *<sup>g</sup><sup>σ</sup>* and *ε*◦ *<sup>f</sup> <sup>σ</sup>* root mean square values for the same *dna* and *dμ* parameters is almost the same: the *ε*◦ *<sup>g</sup><sup>σ</sup>* are also ≈10–30% lower than the *ε*◦ *<sup>f</sup> <sup>σ</sup>*. This means that, for the recognition of decreases small in duration and magnitude, the use of a quasi-Gaussian filter is more preferable than the use of a common FIR filter.

**Figure 4.** Results of calculating the coefficients of relative errors *ε*◦ *gm*, *ε*◦ *f m*.

*3.2. Testing the Method and the Quasi–Gaussian Filter on Experimental Normalized Poisson Data from the URAGAN Hodoscope*

Testing in this section consisted of determining the performance and capabilities of the proposed method and the quasi-Gaussian filter for recognizing small in duration and magnitude decreases in time intervals for the experimental hourly normalized Poisson data registered by the URAGAN hodoscope, taken from [12].

For analysis, a time interval was selected from 09/02/2017, 20:00 UTC to 09/18/2017, 15:00 UTC, with a total duration of 15.6 days. During this interval, the heliosphere was turbulent due to strong solar coronal mass ejections (CMEs) The CMEs that occurred on that period, caused intense geomagnetic storms that were discussed, for example, in [19,20]. The emerging CMEs caused modulations of muon fluxes recorded in MH and led to lower mathematical expectations (including the ones due to Forbush decreases) in Poisson MH data.

MH data were the matrix series of distribution functions of the intensities of muon fluxes *YE*(*i*, *j*, *T*0*n*), defined in a rectangular region *i* = 1, ... ,*N*1, *j* = 1, ... ,*N*2, *N*<sup>1</sup> = 90, *N*<sup>2</sup> = 76, *n* = 1, 2, .... Solid angles correspond to azimuth and zenith indices *i*, *j*, *ϕ<sup>i</sup>* = Δ*ϕ*(*i* − 1), *ϑ<sup>j</sup>* = Δ*ϑ*(*j*−1), Δ*ϕ* = 1◦, Δ*ϑ* = 4◦ in which the registered particles were counted. MH data *YE*(*j*0, *i*0, *T*0*n*) were a time series with indices *j*0, *i*0; the considered interval was determined for *nE* min ≤ *n* ≤ *nE* max, *nE* min = 5900, *nE* max = 6275 (counting hours for [12] began from the first hour of 2017).

Figure 5 shows the results of quasi-Gaussian filtering and interval recognition with reductions in mathematical expectation. The original data *YE*(*T*0*n*) for *j*<sup>0</sup> = 30, *i*<sup>0</sup> = 31 were denoted by light gray jagged lines (index 1). Fluctuations in data with a period of ≈24 h and an amplitude of ≈0.0037–0.0040 are due to the daily rotation of the MH with the Earth. Line with index 2 depicts the data *XEG*(*T*0*n*) filtered based on quasi-Gaussian filter. The recognized intervals of intensity decrease, intensity recovery and intensity mathematical expectation decrease were denoted by a piecewise linear spline-like dashed line *XES*(*T*0*n*) (index 3). Analysis of intervals 5969–6043, 6127–6189 based on *XES*(*T*0*n*) leads to a conclusion that the mathematical expectation values of decreases on them were Δ*m*<sup>1</sup> = 0.01, Δ*m*<sup>2</sup> = 0.005 for the relative decrease rates *dμ*<sup>1</sup> = 0.027, *dμ*<sup>2</sup> = 0.020. For *YE*(*T*0*n*) and *XEG*(*T*0*n*), the mathematical expectations on these intervals were *m*◦ *<sup>E</sup>*<sup>1</sup> = 0.3505, *m*◦ *<sup>E</sup>*<sup>1</sup> = 0.3510, and *m*◦ *EG*<sup>1</sup> = 0.490, *m*◦ *EG*<sup>2</sup> = 0.3520, respectively, on average. The errors of the mathematical expectations estimates were Δ*m*◦ = 0.0010 − −0.0015, which is 10–30% from the mathematical expectation values obtained, and this led to successful recognition of decreases with the relative decrease rates of 0.02–0.03.

Testing on experimental MH data made it possible to draw a conclusion about the efficiency of the quasi-Gaussian filtering method and its satisfactory capabilities for recognizing small fluctuations of the mathematical expectations.

**Figure 5.** Results of quasi-Gaussian filtering and identification of regions with Forbush decreases: original data (line 1), filtered data (line 2), recognized intervals of various muon flux intensity (line 3).

#### **4. Discussion**

The comparison between the model data filtering result obtained using the proposed filter and the one obtained using the fir1 (plots on Figure 3) shows that the resulting time series are close to each other; however, the *XMG*(*T*0*i*) seems to be closer to the initial model. The main quantitative result of testing the method and the quasi-Gaussian filter on model normalized Poisson datasets included the calculations for (14) for a set of realizations/ The resulting errors *σMG* for *XMG*(*T*0*i*), on average, by 15–30% less errors *σMF* for *XMF*(*T*0*i*). This means that the proposed filtering method provided better filtering (noise reduction) than the standard FIR filter. Consideration of the results of statistical modeling made it possible to draw a conclusion about the efficiency of the method and the quasi-Gaussian filter.

Further tests of the new method on model data, aimed at estimating the mathematical expectation and its RMS errors with respect to the durations and magnitudes of model decreases, showed the method capabilities in disturbance recognitions. It can be seen on Figure 4 that the coefficients *ε*◦ *gm* turned out to be less than the values of the coefficients *ε*◦ *f m*, on average, by about 10–30%. The nature of the plots of coefficients *ε*◦ *<sup>g</sup><sup>σ</sup>* and *ε*◦ *<sup>f</sup> <sup>σ</sup>* for the RMS for the same parameters *dna*, *dμ* is almost the same—the coefficients *ε*◦ *<sup>g</sup><sup>σ</sup>* are less than the values of the coefficients *ε*◦ *<sup>f</sup> <sup>σ</sup>*, on average, also by ≈10–30%. From the point of view of recognizing decreases in duration and magnitude, the use of a quasi-Gaussian filter is more preferable than a common FIR filter.

Finally, tests made on real experimental datasets from a muon hodoscope display the method application to data processing and recognition of intervals of decreasing and recovering muon flux intensity. Due to the noise reduction in *XEG*(*T*0*n*), it became possible to clearly see the intervals of quiet data ( Figure 5), intervals with decreases and recoveries and intervals with declines in mathematical expectation; all these recognized intervals were denoted by a line *XES*(*T*0*n*) (index 3 on Figure 5). On the intervals with the boundary points 5900–5954, 6057–6121, 6197–6276 there were quiet data, on the time intervals 5969–6043, 6127–6189 a decrease in mathematical expectation was observed, the time intervals 5955–5970, 6044–6056, 6122–6126, 6190–6196 corresponded to data with decreases and recoveries. On the intervals 5969–6043, 6127–6189, it is quite possible to recognize relative reductions in mathematical expectation. The errors of the mathematical expectations estimates were Δ*m*◦ = 0.0010–0.0015, which is 10–30% from the mathematical expectation values obtained, and this led to successful recognition of decreases with the relative decrease rates of 0.02–0.03 and an average duration of *mathrm approx* 10 h.

Testing the proposed method and quasi-Gaussian filter for data variants with indices *j*<sup>0</sup> = 31, *i*<sup>0</sup> = 30, allowed to obtain results that are almost similar to those depicted on Figure 5); the errors in the estimation of the boundary points of the sections during recognition with depressions amounted to *δn* ≈ 2–5 h. Thus, testing on experimental MH data allowed us to make a conclusion about the efficiency of the method and the quasi-Gaussian filter and their satisfactory capabilities for recognizing mathematical expectation small in duration and magnitude.

#### **5. Conclusions**

The proposed filtering method for time series of normalized Poisson-distributed data, which was based on the developed digital low-pass quasi-Gaussian filter with a finite impulse response, a gain equal to one at low frequencies and non-negative weighting coefficients, turned out to be efficient; the FR of the low-frequency quasi-Gaussian filter of small dimension was characterized by a better approximation to the prototype filter FR than the FR of common FIR filters.

Testing the filtering method based on the quasi-Gaussian filter for the problems of recognizing small in duration and magnitude fluctuation decreases (increases) in mathematical expectations using statistical modeling and statistical tests have confirmed its effectiveness:


Testing the method and the low-frequency quasi-Gaussian filter on experimental Poisson data made it possible to draw a conclusion about its satisfactory capabilities for recognizing decreases with relative decrease coefficients ≈0.020–0.030.

The proposed method of noise reduction and a quasi-Gaussian filter have favorable prospects of using radiation particle counters for digital information processing in problems of experimental physics and measuring technology.

**Author Contributions:** Conceptualization, V.G., A.S. and I.Y.; methodology, V.G.; software, V.G. and V.C.; validation, M.D., A.G. and A.S.; formal analysis, V.C.; investigation, A.S.; resources, A.G.; data curation, A.D., A.K. and N.O.; writing—original draft preparation, V.G.; writing—review and editing, R.S. and M.D.; visualization, V.G.; supervision, I.Y.; project administration, A.S.; funding acquisition, A.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the Russian Science Foundation (project No. 17-17-01215).

**Data Availability Statement:** Data sharing is not applicable to this article.

**Acknowledgments:** The results of experiments presented in this research rely on data collected by the Scientific & Educational Centre NEVOD, National Research Nuclear University MEPhI. We acknowledge URAGAN muon hodoscope data provided by the NEVOD institution. This work employed facilities and data provided by the Shared Research Facility "Analytical Geomagnetic Data Center" of the Geophysical Center of RAS (http://ckp.gcras.ru/) (accessed on 19 April 2021). We would like to thank two anonymous reviewers whose valuable comments helped to improve the manuscript and properly demonstrate the results of our research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

