Real time analysis of the ion density measured by the satellite DEMETER in relation with the seismic activity

. This paper is related to the study of the ion density recorded by the low altitude satellite DEMETER. In a ﬁrst time there is an automatic search for ionospheric perturbations in the complete satellite data set of ion densities. Then perturbations due to known ionospheric phenomena (for example, solar activity) are eliminated as well as perturbations not above a seismic zone. In a second time, there is a search to know if each selected perturbation corresponds to a future earthquake. The earthquakes have been classiﬁed depending on their magnitude and depth. This attempt to predict earthquakes of course generates false alarms and wrong detections. The results of this statistical analysis are presented as function of various parameters. It is shown that the number of false alarms is very important, because the ionosphere has variations not only linked to the seismic activity.


Introduction
Ionospheric disturbances can be considered as short-time precursors. although the appearance of changes in the ionosphere is not firmly established before all earthquakes (Liperovsky et al., 1992;Molchanov, 1993;Hayakawa et al., 1996;Hayakawa, 1997Hayakawa, , 2007Liu et al., 2000;Chuo et al., 2002;Singh et al., 2003;Akhoondzadeh et al., 2010;Bankov et al., 2010;Kon et al., 2011). As early as in May 1960, during a period of quiet geomagnetic activity, the critical frequency of the F2 layer (foF2) was slightly below its average value a few hours before the largest Chilean earthquake recorded in history which occurred on 22 May 1960 (Foppiano et al., 2008). Irregularities in the ionospheric sounder data were also found before the Alaskan earthquake taking place on 28 March 1964 (Davies and Baker, 1965).
As a result of researches during the last ten years, it has been shown that the ionosphere is unexpectedly extremely sensitive to the seismic effect (Hayakawa and Molchanov, 2002). There have been recently accumulated reports on the ionospheric precursory effects of earthquakes (EQs) especially after the launch of the DEMETER satellite in 2004. The abnormities in TEC (Total Electron Content) were observed 4-6 days before the L'Aquila 6.2 earthquake with a maximum increase of about 67 % (Dabas et al., 2007;Akhoondzadeh et al., 2010;Stangl et al., 2011). Hasbi et al. (2011) investigated the ionospheric variations before four earthquakes with M ≥ 7.9 that occurred during [2004][2005][2006][2007] in Sumatra using GPS and CHAMP data, and the result shows the occurrence of positive and negative anomalies detected within a few hours to 6 days before the earthquakes. Hayakawa et al. (2008Hayakawa et al. ( , 2011 analyzed subionospheric VLF propagation data and found perturbations before the 2007 Niigata Chuetsu-Oki 6.8 earthquake and the 2010 Haiti 7.0 earthquake with an abnormal radius up to several kilometers. On 9 May, just 3 days before the Wenchuan 8.0 earthquake,

Published by Copernicus Publications on behalf of the European Geosciences Union.
TEC and foF2 of GPS observations, electron density and ion density of DEMETER satellite were subjected to an abrupt increase or decrease with amplitude being up to 40 %-80 %. The spatial distributions of the anomalous and extreme reductions and enhancements indicate that the earthquake preparation area is about 1650 km and 2850 km from the epicenter in the latitudinal and longitudinal directions, respectively (Zhao et al., 2008;Yu et al., 2009;Liu et al., 2009;Zeng et al., 2009;Lin et al., 2009;Akhoondzadeh et al., 2010;Zhang et al., 2009;Xu et al., 2010Xu et al., , 2011. Seismo-ionospheric disturbances within a few days before earthquakes have been registered in 73 % of earthquakes with magnitude 5.0, and in 100 % of earthquakes with magnitude 6.0 (Pulinets et al., 2003). Silina et al. (2001) analyzed Es-spread effect before 36 earthquakes with magnitude 5.5 ≤ M ≤ 6.6, which are classified into "deep" (h > 33 km) and "crust"(h ≤ 33 km) groups. The result shows that Esspread increases 1-3 days before earthquakes, the effect being stronger for deep ones. Using TEC data, Liu et al. (2004) statistically described the temporal parameters of the seismoionospheric precursors detected 1-5 days prior to 20 major earthquakes in Taiwan. Akhoondzadeh et al. (2010) statistically analyzed the variations of the electron and ion densities prior to four large earthquakes using IAP (Instrument Analyseur de Plasma) and ISL (Instrument Sonde de Langmuir) experiments onboard the DEMETER satellite and with GPS measurements. The appearance of positive (increase of the density) and negative (decrease of the density) anomalies happens during 1 to 5 days before all studied earthquakes during quiet geomagnetic conditions. Parrot (2011) statistically analyzed the ion density perturbations of the DEME-TER satellite related to the seismic activity. Starting from an earthquake database, an automatic search for perturbations has been done on data recorded by the satellite close to the epicenters until 15 days before the occurrence of the earthquakes. The result shows that the number and the intensity of the ionospheric perturbations are larger prior to earthquakes than prior to random events, and that the perturbations increase with the magnitude.
In this paper a statistical analysis with the whole ion density data recorded during about 6.5 yr will be shown. The DEMETER payload is briefly described in Sect. 2. In Sect. 3, the data processing method is presented. It uses an automatic detection of peaks in the ion density recorded along orbits and a search to know if peaks correspond or not to earthquakes. The results of the statistical analysis are presented in Sect. 4. Discussion and conclusions are provided in Sect. 5.

The DEMETER satellite
DEMETER is a low-altitude satellite (710 km) launched in June 2004 onto a polar and circular orbit which measures electromagnetic waves and plasma parameters all around the globe except in the auroral zones (Parrot, 2006). The altitude of the satellite decreased to 660 km in December 2005. The satellite's science mission came to an end in December 2010. Due to technical reasons (onboard attitude control system which generated interferences), data were only recorded at invariant latitudes less than 65 • . The orbit of DEMETER is nearly sun-synchronous, and the up-going half-orbits correspond to nighttime (22:30 LT), whereas the down-going half-orbits correspond to daytime (10:30 LT). Nearly sunsynchronous means that, everyday, the satellite does not return exactly above the same point but above the same area (it could be at more than 1000 km from the point flown over the days before). To study the variations of the ion density, the payload of DEMETER includes the instrument named IAP which gives ion density with a 4-s time resolution. Details of the IAP experiment can be found in Berthelier et al. (2006).

Description of the data processing
In a way the current data processing of this work is complementary to the method used in Parrot (2011) where ionospheric perturbations were automatically searched in the vicinity (time and space) of known earthquakes. Here all ionospheric perturbations are automatically searched in the complete data set, and then it is automatically checked if a perturbation could correspond to a given earthquake. As the DEMETER micro-satellite runs to the end in December 2010, the complete data set for the parameter IAP is used. In order to process this amount of data step by step, we have designed two software. The aim of the first software is to search for ion density (the sum of H + , He + and O + ) perturbations everywhere around the Earth. On one hand, a map of seismic activity zones is constructed using the list of earthquakes with M W ≥ 4.8 taking place from 1 July 2004 to 31 December 2010, which basically covers the whole period of the satellite's life. On the other hand, a complete set of half-orbit IAP ion density files is used as the original input. There are burst mode and survey mode in the raw data, and the sample rate of the former is two times the latter. The average of two adjacent points for the burst mode is calculated in order to make the input data with the same sample rate. The SAVGOL method is employed to smooth the data before searching for perturbations. The SAVGOL function returns the coefficients of a Savitzky-Golay smoothing filter (Savitzky and Golay, 1964). This filter is defined as a weighted moving average with weighting given as a polynomial of a certain degree, and it can use any number of points for this weighted average. A detected perturbation is kept if it meets two basic limit conditions: 1. It must have a definite duration time. For example, we set the maximum time at 2 min (this was tested many times and considered to be reasonable for most events; taking into account the satellite velocity, this corresponds to 840 km), and the minimum time at 23 s (about 5 data points) to avoid spurious impulsions caused by one or two points.
2. If the peak (decrease or increase) does not lie in or near a seismic zone (defined by the earthquake map described above), it should be eliminated. For instance, in the course of our following data processing, if the distance between a perturbation and its nearest seismic zone is larger than 1500 km, this peak will not be kept as it is considered to have no relation with any earthquakes.
All these parameters can be manually changed in the software. Values of the peaks are compared to background values. For example in case of increase, a maximum is first determined. This maximum being at the time t, a search is done for the two minima that are just around t (first change of the sign of the derivative on each side of the maximum). If they occur at t1 (< t) and t2 (> t), the perturbation is kept if 23 s < t2-t1 < 120 s, and the background value is determined by values at t1 and t2.
At the end, all the accepted perturbation information, such as day, time, orbit, sub-orbit, background value, amplitude, trend (increase or decrease, which means the trend when the amplitude is compared with the background value) and the peak's last time, forms an output database file. The location of every perturbation can be also drawn in a global map (Fig. 1). The aim of the second software is to check whether the events are corresponding to an earthquake or not. The interface is shown in Fig. 1. In this step, the output database file from the first software is used as input file. The K p index is taken into account in order to reduce the effect of the geomagnetic activity on the ionosphere. This geomagnetic activity is induced by solar magnetic storms. The input parameters also include the minimum magnitude M W , the distance (D) between the location of the perturbation on the orbit and the epicentre, the depth of the seismic epicentre below the Earth's surface (d), and the delay time before and after an . For each perturbation of the list, we need to check if this perturbation could correspond or not to one earthquake under the limit conditions above. If an earthquake is corresponding to one or to more than one perturbation, we consider it is a good detection; if not, it is a bad detection. If a perturbation corresponds to an earthquake, it is a right alarm; if not, it is a false alarm. These numbers form the analysis report file at the end. The locations of the earthquakes for the good detection can be shown on the global map. In addition, the software can automatically perform the statistics concerning the number of the perturbations corresponding to the same earthquake (good detection) and the trend ("Increase" if the maximum value of the peak is larger than the background value; "decrease" on the contrary) and then output the result. If any one of the limit parameters is changed, we will get different results.

The statistical analysis of the results
According to the results of the first software, the number of the total IAP ion density data files from 11 August 2004 to 16 December 2010 is 96 863, including 27 257 933 data points with the same sample rate. The number of detected perturbations is 56 139 with D being 1500 km (in fact here D is equivalent to the maximum distance between the location of the perturbation and a seismic area) and the maximum time being 2 min. The earthquake database considered in the present paper includes 21 863 earthquakes from 20 August 2004 to 31 December 2010, with M W ≥ 4.8 (USGS: http://www.usgs.gov). To check if the earthquake magnitudes have an influence, they are classified into three groups according to the magnitude M W (see Table 1). Table 1 shows that the number of earthquakes with M W = 4.8-5.0 is 12 057, 55.1 % being more than half of the total earthquakes, 41.0 % with M W = 5.1-6.0, and 3.9 % with M W larger than or equal to 6.1.
Not all ionospheric perturbations are affected by earthquakes. They may be due to other sources, such as solar activity, acoustic gravity waves (AGW), travelling ionospheric disturbances, plasma dynamics and large meteorological phenomena (Parrot, 2011). We keep K p < 3 in the course of the following data processing to try to eliminate perturbations due to solar activity. The possible effects of the other phenomena are smoothed when we consider a huge number of events.
According to the definition of all parameters, we set different ranges for each of them. First, earthquakes are divided into two groups on the base of the epicentre depth d. All earthquakes with d ≤ 20 km are selected to be one group in order to statistically study the ionospheric influence of the epicentre depth. For the group of 4.8 ≤ M W ≤ 5.0, the number of earthquakes with d ≤ 20 km is 4487 accounting for 37.2 %; 3783 and 42.3 % for 5.1 ≤ M W ≤ 6.0; 403 and 47.2 % for M W ≥ 6.1. Then, for each group of earthquakes, relations between perturbations and earthquakes are searched according to different D values (500 km, 1000 km and 1500 km) to analyse the effects of the distance. Considering the time of a given peak, we also search up to 7 Nat. Hazards Earth Syst. Sci., 12, 2957-2963, 2012 www.nat-hazards-earth-syst-sci.net/12/2957/2012/ and 15 days after for the finding of an earthquake at a distance < D.
For each group of earthquakes, the result changes when the limit conditions D, d and T vary. The complete results for the parameters, D = 0-500 km, d = 0-20 km, T = 0-7 days, and M W = 4.8-5 are the following: Number of perturbations meeting the limit conditions = 46 446; number of earthquakes meeting the limit conditions = 4487; number of right alarms = 767; number of false alarms = 45 679; number of bad detections = 3880; number of good detections = 607. The number of right alarms does not match the number of good detections, because, for a given earthquake, we can have several perturbations. Concerning the number of false alarms, it must take into account the number of right alarms for earthquakes with magnitudes larger than 5 in the same conditions (727 for 5.1 ≤ M W ≤ 6.0 and 89 for M W > 6.1). It means that the real number of false alarms is 44 863. The results for several values of the input parameters (including the previous ones) are shown in Tables 2-7. In these results, n designates the average number of perturbations for each earthquake detected, r the ratio of the number of earthquakes with good detections and the number of earthquakes that comply with the limit conditions, and s is the ratio of the number of right alarms with the trend being "increase" and the total number of right alarms.

Discussions and conclusions
Checking the results in the tables, one can see the following: -When we increase the time to search for earthquakes after the occurrence of perturbations from 0-7 days to 0-15 days, the ratio r increases which is obvious as we increase the possibilities to find an earthquake.
-When we increase the parameter D from 0-500 km to 0-1000 km (or to 0-1500 km), the ratio r increases for the same obvious reasons.
-Most of the perturbations are increases.
Overall the results are far from being good because of the following reasons: -The number of false alarms that are due to natural variations of the ionosphere is always very high. In a future work, we will try to select perturbations with special parameters (duration, amplitude relative to the background) in order to reduce this number.
-The number of wrong detections is also high. This could be due to the fact that we use a satellite that does not allow continuous measurements above a given seismic zone. We do not expect to have continuous ionospheric perturbations, and the satellite is "above" a seismic area only a few minutes per day. Then, we can miss a perturbation. It is possible to decrease this number of wrong detections if several satellites are simultaneously used.
However, it remains that a relation between the perturbations and the earthquakes exists because of the following: -The ratio r between the number of detected earthquakes and the total number of earthquakes increases with the magnitude.
-When we increase the depth of the earthquakes, this ratio r decreases.
These last two points are intuitively expected. At the end, one must say that a good detection does not mean that we are able to predict the corresponding earthquake. Uncertainty about its position is large. There is a way to find more precisely its position if we have several perturbations for one earthquake. It is possible to use triangulation and also to use the position of the fault below the satellite. But uncertainties about the time and the magnitude are also important.