2.1 Types of Particle Detector Failures Encounter during Cosmic Ray Flux Multiyear Monitoring

2.2 Moving Median Filter (MMF) – “Horizontal” Median

2.3 Relational Median Filter for Multichannel Measurements (RMF); “Vertical Median”

2.4 Verification of Filtering Algorithms

2.5 Monitoring of the Stability of Measuring Channels

2.6 Conclusion

The influence of the Sun on GCR flux can be described as a modulation of the usually stable “background”. The Sun modulates GCR in several ways. The explosive flaring processes on the Sun result in the ejection of huge amounts of solar plasma and in acceleration of the copious electrons and ions. These particles comprise, so called, Solar Cosmic Rays (SCR). The SCR reach the earth and initiate secondary elementary particles in the terrestrial atmosphere, increasing the count rates of particle monitors by several percents. This effect is called Ground Level Enhancement (GLE). The solar wind “blows out” the lowest energy GCR from the solar system, thus changing the GCR flux intensity inverse proportionally to the Sun activity. Huge clouds of magnetized plasma, so called Interplanetary Coronal Mass Ejections (ICMSs) usually proceeded by shocks are travelling in the interplanetary space tens of hours before reaching earth. During their travel ICMEs interact with GCRs changing their intensity and direction. Particle detectors on Earth register these “modulation” effects of ICME as changes of intensity of the secondary cosmic rays. On arrival at the earth the magnetic field of the plasma cloud deplete the GCRs, measured as decrease in intensity of the secondary cosmic particles (so called Forbush decrease).
Worldwide networks of particle detectors are continuously monitoring changing fluxes of the secondary cosmic ray species. Neutron monitors (NM, Moraal et al, 2000) are measuring neutron fluxes, Solar Neutron Telescopes network (Tsuchiya et al., 2001) is searching very rare neutron fluxes from the Sun, Muon network (Munakata et al., 2000) is measuring high energy muon flux. The SEVAN network (Chilingarian et al., 2009) is measuring neutral particle fluxes and, additionally, charged particle fluxes with different energy thresholds. The data from networks is collected in special data bases (NMDB, FP7) and usually is available on-line. The, so called, time series of the count rates of particle detectors comprise the detailed (not worse than 1-minute) data file of the upmost important geophysical information for the research of solar-terrestrial connections and for issuing forewarnings and alerts on the violent conditions of the Space Weather.
Particle detectors, usually, located in difficult of access regions sometimes malfunction due to altering of hardware during long-year exploitation and it is not easy to perform fast repair. Therefore the correction of raw data (filtering, smoothing) is necessary step of data analysis. Filtering algorithms are usually based on the comparisons of data from identical measuring channels. Particle detectors of the world-wide network of Neutron Monitors (Clem & Dorman et al., 2000) usually consist of 3 sections, 6 identical proportional counters in each. If the ratio of count rates of different sections is changing within defined limits the detector overall count rate is performed by simple summon of all sections.

2.1 Types of Particle Detector Failures Encounter during Cosmic Ray Flux Multiyear Monitoring

Mainly there are 4 kinds of errors, as shown in Figure 1-4and different combinations of these types.

Figure 1. Abrupt Spike (jump).	Figure 2. Slow Drift.

Figure 3. Abrupt change of mean continued with recovery.	Figure 4. Abrupt change of mean without recovery.

Mentioned errors in particle detector’s operation can lead to erroneous estimation of the Fd magnitude; prevent detection of the GLE, (usually not very large at middle latitudes - 1.5-2%), etc…. We introduce algorithms based on the stabilizing properties of the median for correction of multichannel detector’s data. We heavily use the “overabundance” of ASEC data due to numerous identical channels measuring one and the same physical quantity (flux of particles of definite type).

2.2 Moving Median Filter(MMF)

– “horizontal” median

Window of size L is “moving” from beginning to the end of time series with unit step; the median of L values falling within window is continuously calculated. If difference of the current time series and median is too large several actions are triggered:

a) the outlier can be substituted by the median value;

b) the outlier can be substituted by the special code;

c) execution of program stops and an request to operator is send.

Also MMF performs time series smoothing:

d) All time series are substituted by the corresponding moving median values.

If number of outliers exceed pre-chosen limit another filtering algorithm described in the next section is invoked.

In Figure 5 the schematic chart of MMF is presented.

Figure 5. Schematic view of the “horizontal” window "moving" across time series.

Algorithm Description:
Notion:

• Time series of detector channel at moment i is denoted by small letter vi; median of the L successive elements of time series started at i is denoted by
• Moving window width – L;
• Minimal and Maximal values of the wsindow width – L_minand L_max L_min < L_max;
• Maximal and minimal possible value of time series median – P_max,P_min;
• Maximal possible deviation of time series from median value D_max.

Algorithm steps:

1. Select time series from database with N elements;
2. Start algorithm operation from the first of time series, assign i = 0;
3. Define L=L_min; if i < N, then assign i=i+1 and continue; otherwise write filtered time series into data base; calculate length of periods when algorithm substitute the time series by the median value, send to operator all messages and stop.
4. Select L-1elements of time series to the right; calculate the median value M_i,Lif its value is in the limit of the predetermined values M_i,Lε(Pmin –Pmax ) then continue; otherwise, check if L < L_max enlarge L by 2 and repeat steps 3,4; otherwise report about algorithm failure at point i and store algorithm parameters for ith time series: (i, V_i, L_max); then go to 3
5. Check if abs(V_i - M_i,L)<D_max then continue; otherwise, report erroneous i-th time series, store algorithm parameters (i,V_i,M_i,L) and assign V_i=M_i,L then go to 3.

We describe algorithm operation for the option a), options b),c),d)) can be readily obtained by the minor changes of the described algorithm.
This filter is optimal for the time series containing spurious short spikes and abrupt changes of mean followed by recovery Figure 2 and Figure 4. It cannot correct smooth trends due to slowly changing parameters of particle detectors due to altering of Photomultiplier (PM) or/and electronics elements of Amplitude-to-Digital-Converter (ADC) and power supply (Figure 3, 5). Using MMF we can smooth time series, and get rid of short spikes, but it can’t solve the problem when time series mean is changing gradually during large time span; another algorithm should be used to correct such defects using the data of other channels of the same monitor.

2.3 Relational Median Filter for Multichannel

Measurements (RMF); “Vertical Median”

Let's suppose that detector consists of M identical channels, however due to individual characteristics of sensors used (photomultipliers, proportional counters, etc…) the mean count rates of channels , j=1, M are dispersed within definite (not very large) limits.
Notion:
M– number of channels of the monitor;
- mean count rate of j-th channel;
N_total – sum of mean values of all channels (detector mean count rate);s
– Median value of M channels at i-th minute;
F_j-the equalizing coefficient of j-th channel;
V_i^j - i-th time series of j-th channel; V_i – estimate of the total detector count rate at moment i.
At the start of detector operation by assigning to each channel the appropriate coefficient F_j; j=1, M it is possible to equalize the mean count rates:
(1)
Only after this “equalizing” operation, it is valid to calculate the median. The detector count rate at moment i can be calculated according to:
(2)
The median estimate of count rate is much more stable in presence of outliers (bad channels) though its variance is greater comparing with mean value in absence of outliers.
Also if jth channel of detector is continuously and incoherently changing (operating unstable according to reports of MMF) its time series can be substituted by the median value:

(3)

The possible scenario of implementation of both algorithms can be as follows:<\p>

1. For some initial period of detector operation possibly without any errors the mean count rates and coefficients F_j are calculated and stored.
2. At the end of the day the data of all channels of detector are filtered with MMF algorithm;
3. If some channels operate unstable according to reports of the first algorithm RMF turns on, it reads the stored means and coefficients and corrects the malfunctioned channel data.
4. Channel means and appropriate coefficients are renewed and stored.
5. If second algorithm did not correct the data (which means that all or nearly all channels have been corrupted or detector was switched off due to some overall failure) system sends an e-mail to manager.

By automatically implementing both filtering algorithms each day and storing renewed mean values of channels and appropriate coefficients it is possible to correct all mentioned failures (demonstrated in Figure 2,3,4,5). The time history of the equalizing coefficients will help to outline non-stable channels and repair them. Below we demonstrate some examples of filtering by these algorithms.
In Figure 7 we present the corrected data of Aragats Neutron Monitor in May 2008 (see Figure 1 for raw data). Time series in the Figure are artificially shifted from each other for better assessment. The combination of MMF and RMF corrects data, eliminating spikes and filling corrupted periods. The Data in all graphs shown here has been smoothed by constant window L=60. The data of Nor-Amberd Neutron Monitor during 23rd solar cycle was smoothed by changing width of the window started from L=60 to L_max=600 (if needed).

Figure 6. Aragats Neutron Monitor, MAY 2008 – after filtering.

In Figure 12 we present the Nor Amberd neutron monitor data as measured during 23rd solar cycle (1997-2007). Due to numerous failures of several detectors in the end of cycle, there are significant biases in the count rate of monitor. Efficiency of several channels due to failures of high voltage supplies and aging effects of counters itself go down and overall count rate of monitor also go down. As we can see in Figure 18 after applying filtering algorithms overall pattern of monitor count rate changes according the solar activity and comes very close to data of Alma-Ati Neutron Monitor located at altitude and latitude close to Aragats ones.

Figure 7. The pressure sensor time series at Nor-Amberd Station, before filtering .	Figure 8. Pressure at Nor-Amberd Station, after filtering.

In Figure 10, Figure 11, we present correction algorithm operating according to equation (2). After finishing of its execution, the time series the MMF algorithm started operation. Although the variance of the filtered time series is larger comparing with initial ones, it corrects many mistakes apparent in the raw data (black – raw data, magenta – corrected).

Figure 9. Correction of the AMMM time series,>5 Gev Muons.

Figure 10. Correction of the ARNM time series.

In Figure 12 we present the Nor Amberd neutron monitor data as measured during 23rd solar cycle (1997-2007). Due to numerous failures of several detectors in the end of cycle, there are significant biases in the count rate of monitor. Efficiency of several channels due to failures of high voltage supplies and aging effects of counters itself go down and overall count rate of monitor also go down.
As we can see in Figure 18 after applying filtering algorithms overall pattern of monitor count rate changes according the solar activity and comes very close to data of Alma-Ati Neutron Monitor located at altitude and latitude close to Aragats ones.

Figure 11. Cosmic Ray intensity modulation by NANM, 23rd solar cycle – before filtering.

<

Figure 12. NANM, 23rd solar cycle – after filtering.

2.4 Verification of Filtering Algorithms

To check the suggested filtering techniques we perform simulation study with artificially corrupted time series. 18 time series with 1 million points have been simulated according to mean values and variances of the channels of Nor-Amberd Neutron Monitor. A trend has been introduced to time series to imitate solar modulation. Then all 3 described types of failures were introduced in the time series (see Figure 15). In Figure 16 one can see that after correction with median algorithms all errors have been washed out and the solar “modulation” effect is apparently seen.

Figure 13. Simulated Time Series

Figure 14. One of 18 simulated time series with modulation effect – before filtering.

Figure 15.One of 18 simulated time series with modulation effect – after filtering.

2.5 Monitoring of the Stability of Measuring Channels

During multiyear operation of particle monitors, mean count rates continuously alter not only by solar modulation or possible entering of regions where Galaxy arms are sending abundant GCR from supernovae explosions, but also by such prosaic effects as electronic components aging. Therefore, to identify instrumental failures and to avoid exploration of the artifacts instead of new physics we have to monitor carefully and continuously detector parameters. It is possible to do it by monitoring equalizing coefficients of monitor channels. The monthly (or decade) plots of coefficients will help to find unstable channels.
The channel mean count rates are changing due to solar modulation effects, in contrast, the equalizing coefficients (see section 3) should be stable despite changing means.
Therefore, it will be much easier to detect non-stable channels by monitoring the plots of coefficients, than to change channel means. Figure 17and 18 are an example of our approach.

Figure 16. Day-to-day changes of the mean values of Aragats Neutron Monitor; at November 22 there were power supply cut off.

Figure 17. Day-today changes of the channel coefficients of Aragats Neutron Monitor.

Although from Figure 17 we can notice that the variations of two channels are significantly larger than variations of the 16, in Figure 18 the behavior of the corresponding coefficients demonstrate failure of two channels much more pronounced. The same method can be implemented to check data from numerous detectors. Several Eurasian countries joint the efforts to establish Neutron Monitor Data Base (NMDB, data available from the NMDB.eu). Data from numerous neutron monitors is gathered in NMDB and for the “housekeeping”, it is necessary to run periodically tasks for checking if all monitors are adequate and data is correct and coherent. We select time series from several monitors participating in the NMDB project to run the “coherence” test (see Table 1).

Table 1. Neutron Monitors used in the “coherence” test

Neutron Monitor	Altitude,m	Rigidity,GV	Type
Alma Ata	3340	6.69	18NM64
Rome	60	6.32	17NM64
Aragats	3200	7.14	18NM64
Nor-Amberd	2000	7.14	18NM64
Moscow	200	2.46	24NM64
Oulu	0	0.81	9NM64
Athens	260	8.53	6NM64

We have taken pressure corrected data of these monitors for time period 24.10.2008-25.12.2008 and calculated coefficients for these seven time series, according to equations 1-3 . In Figure 19 the raw data are posted; in Figure 20 we present the median corrected, according to equation (3), data of NMDB Monitors (note, that the spikes in Nor Amberd and Izmiran monitors are filtered out and gaps are filled). In Figure 21 we present the “equalizing” coefficients for all seven monitors calculated for each day of selected period. As you can see, 6 coefficients from 7 demonstrate very stable behavior, proving that all parameters of the neutron monitor chambers remain stable and constant. The calculated coefficients for the Athens monitor are much more variable. The high variability (non-coherence) of the Athens monitor’s coefficient may be caused by drift of the electronics parameters (including pressure sensor) in the end of 2008. When dealing with multiple remote sensors it is of vital importance to develop a number of quality tests to check continuously the data coming from different remote destinations. Although data of the NMDB network detectors are similar, different groups use different data acquisition electronics, pressure sensors and data transfer protocols. The time history of the “equalizing coefficients” is one of such tests to help keep NMDB data reliable and adequate for further physical inference.

Figure 18. Pressure Corrected Data from NMDB

Figure 19. Pressure Corrected and Median Smoothed Data from NMDB

Figure 20. Equalizing coefficients of the NMDB facilities

2.6 Conclusion

Filtering of the multichannel data of particle detectors having operated for many years for the detection of the solar modulation effects and, maybe, sidereal modulation effects, is of vital importance. During multiyear measurements, characteristics of detector undergo critic changes due to aging effects of sensors and discrete elements of electronics. Overabundance of the information allows introducing correction algorithms using stabilizing properties of the median of time series. Continuous storing and monitoring of the mean values of all channels along with their equalizing coefficient allows archiving of time-history of the behavior of all channels. Examining the relative behavior of channel means and coefficients during multiyear operation it became possible to distinguish the physical effects from instrumentation failures. See for example discussion in (Buetikofer & Flueckiger, 2008) and (Bieber et al., 2007). Our approach also allows not only correction of mistakes due to hardware malfunction, but simple and efficient method of timely detection of non-stable channels or/and mistakes in data bases collecting time series from different remote detectors.

References

• Buetikofer R., Flueckiger E.O., Long-term stability of Swiss neutron monitors. Private communication, 2008.
• Belov A., Blokh Ya., Klepach E., et al. Primary Processing of Cosmic Ray Station Data: Algorithm, Computer Program and Realization. Kosmicheskie Luchi. 25, pp. 113-134, 1988 (in Rusian).
• Bieber J. W., Clem J., Desilets D., Long-term decline of South Pole neutron rates. J. Geophys. Res., 112, A12102, 2007
• Chilingarian A., Hovsepyan G., Arakelyan K., et al., Space environmental viewing and analysis network (SEVAN). Earth, Moon, and Planets, 104, 195, 2009.
• Moraal, H., Belov, A., Clem, J.M.: Design and co-Ordination of Multi-Station International Neutron Monitor Networks, Space Science Reviews 93, 285-303, 2000
• Munakata, K., J.W. Bieber, S. Yasue et al.: Precursors of geomagnetic storms observed by the muon detector network, J. Geophys. Res. 105, 27457-27468, 2000
• Tsuchiya, H., Muraki, Y., Masuda, K., et al. Detection efficiency of a new type of solar neutron detector calibrated by an accelerator neutron beam, NIM, A 463, 183 – 193, 2001.

Attachment	Size
figure2_2.jpg	9.81 KB
figure2_3.jpg	8.02 KB
figure2_4.jpg	7.59 KB
figure2_5.jpg	19.14 KB
figure2_6.jpg	36.89 KB
figure2_7.jpg	50.94 KB
figure2_8.jpg	44.51 KB
figure2_9.jpg	28.45 KB
figure2_10.jpg	31.79 KB
figure2_11.jpg	22.6 KB
figure2_12.jpg	31.77 KB
figure2_13.jpg	50.21 KB
figure2_14.jpg	28.33 KB
figure2_15.jpg	24.34 KB
figure2_16.jpg	34.81 KB
figure2_17.jpg	34.43 KB
figure2_18.jpg	39.9 KB
figure2_19.jpg	46.12 KB
figure2_20.jpg	29.14 KB
2-1.gif	1.63 KB
2-2.gif	831 bytes
2-3.gif	819 bytes
figure2_5.gif	5.12 KB
2-1-0.gif	2.17 KB
figure2-1.jpg	3.43 KB
figure2-3.jpg	8.02 KB
n_j.gif	301 bytes

‹ 1. Modeling of the detector response. up 3. Calculating statistical significance of the detected peaks in time series. ›