A Statistical Analysis of TIR Anomalies extracted by RST in Relation with Earthquake in Sichuan Area with Use of MODIS LST Data

There is a long history for research of earthquake prediction, but weakness of traditional approaches to study seismic hazard have been more and more evident. Remote sensing and earth observation technology, which is a new method that can instantly acquire a large area of abnormal information caused by earthquakes, is believed to be the key to the breakthrough of the bottleneck 10 in the study of earthquake prediction. A multi-parametric approach seems, instead, to be the most promising approach in order to increase reliability and precision of short-term seismic hazard forecast, and Thermal Infrared (TIR) anomaly is an important part of the earthquake precursors. Though many scientists have studied the correlation among TIR anomalies identified by the Robust Satellite Techniques (RST) methodology and single earthquake, there is few study to extract the TIR 15 anomalies in long period and large study area. Moreover, a statistical analysis of TIR anomalies in relation with earthquake is needed to determine whether there is the existence of TIR anomalies before earthquake. In this paper, a refined RST data analysis and Robust Estimator of TIR Anomalies (RETIRA) index were used to extract the TIR anomalies from 2002 to 2018 in Sichuan area with use of Moderate-resolution Imaging Spectro-radiometer (MODIS) Land Surface 20 Temperature (LST), and the earthquake catalog were also used to study the correlation between TIR anomalies and occurrences of earthquake. Most of the thermal infrared anomalies correspond to earthquakes, and statistical methods are used to prove that there is a correlation between the extracted thermal infrared anomalies and earthquakes. And this is the first time to evaluate earthquakes prediction ability with use of PPV, FDR, TPR and FNR, the statistical result shows that 25 the prediction ability of RST in Sichuan area is limited.


Introduction 30
There are numerous observations of surface temperature changes prior to Earth's crust earthquakes (Tronin, Hayakawa et al., 2002).Nowadays, Thermal Infrared Remote Sensing has taken to be a new method for seismic precursors detecting.Anomalous thermal infrared emissions have been widely detected by satellite sensors before the major earthquakes (Piroddi, Ranieri et al., 2014).Several studies discovered space-time anomalies in TIR and outgoing longwave satellite 35 imagery, from weeks to days, before and after earthquakes (Wang, 1984, Gorny, Salman et al., 1988, Qiang, Xu et al., 1991, Tronin, 1996, Tramutoli, Bello et al., 2001, Ouzounov and Freund, 2004, Tramutoli, Corrado et al., 2015).Identification of thermal infrared (TIR) precursors as pre-seismic signal has gained support over world, especially in Russia, China, India, United States, Italy, and Saraf et al. observed such short-term anomalies around the epicentral region for earthquakes in India, 40 passive microwave DMSP-SSM/I satellite data and call them 'transient TIR anomalies (Saraf, Rawat et al., 2009).
There is few data analysis techniques that can isolate residual TIR variations, potentially associated with earthquake occurrence, from the normal variability of TIR signal due to other causes (Tramutoli, Cuomo et al., 2005).But, more than 10 years (since 2001) of application of the 5 general Robust Satellite Techniques (Tramutoli, 1998, Valerio, 2005, Tramutoli, 2007) methodology to this issue, have shown the ability of this approach to discriminate anomalous TIR signals possibly associated with seismic activity from normal fluctuations of Earth's thermal emission related to other causes(e.g., meteorological) independent of the earthquake occurrences (Eleftheriou, Filizzola et al., 2016).RST is based on the Robust AVHRR techniques (RAT), which is proposed for 10 environmental monitoring with the use of NOAA/AVHRR observations (Tramutoli, 1998).And since then most of the announced RAT applications have demonstrated their reliability as well their exportability on different satellite sensors and geographic areas, so RAT evolved into RST (Tramutoli, 2007).RST contains two main steps: the first RST requirement is the characterization of behavior in normal conditions; the second step is the establishment of change 15 detection criteria which should be specified for each considered phenomenon class and chosen technology as well as for the time and place of the observation (Tramutoli, 2007).
With use of RST, many researchers have extracted TIR anomalies (in the following all the TIR anomalies refer to the thermal infrared anomalies extracted by RST) in different earthquakes, and have analyzed the space-time distribution of TIR anomalies.Using MODIS LST data, Pergola 20 studied the 6 April 2009 Abruzzo earthquake found that spatially extended and time persistent TIR anomalies (with RETIRA>3) appear in some space-time correlation with earthquakes of different magnitude occurred in Italy in the considered period (15 March-15 April) and since seven days before the Abruzzo main shock (Pergola, Aliano et al., 2010), while Mebrouk studied 21 May 2003 Boumerdes earthquake and detected a thermal anomaly persisting for a week during the month 25 preceding the earthquake (Bellaoui, Hassini et al., 2017).Many researchers also used data of other satellites, Aliano used 8 years of Meteosat TIR observations to analyze 21 May 2003 Boumerdes/Thenia(Algeria) earthquake found that the area of interest was affected by significant positive thermal anomalies (S/N>2.5-3)about one month before the main shock (Aliano, Corrado et al., 2007), and M.Lisi studied 6 April 2009 Abruzzo earthquake with the use of NOAA/AVHRR 30 TIR observations and TIR anomalies were identified in some space-time correlation with Abruzzo earthquake epicenter between 30 March and 1 April (Lisi, Filizzola et al., 2010).Genzano have studied the 2009 Abruzzo with use of different satellite data(5 years of MSG/SEVIRI, 15years of NOAA/AVHRR and 8 years of EOS/MODIS), t no similar results have been observed in confutation (Genzano, Corrado et al., 2010).Besides the analysis for TIR anomalies in single 35 earthquake, Tramutoli et, al. have studied the causes of TIR anomalies, they performed a test over an area affected by variable gas emissions to study the correlation between TIR anomalies and seismicity, and found that the general gas dispersion models and the spatial features follow the hypothesis of a strict relation between greenhouse gas releases and TIR anomalies related to seismic activity (Tramutoli, Aliano et al., 2013).40 Nowadays, some researchers have done a long-term statistical analysis to determine the correspondence between TIR anomalies and earthquakes.Genzano used data of GMS-5/VISSR TIR measurements to study earthquakes with M>4 occurred in a wide area around Taiwan, in the month with M>4 or M>4.5 were considered, the false positive rate remained less than 6% when the M >5 is applied (Genzano, Filizzola et al., 2015).Tramutoli et al. studied the earthquakes with M>4.0 in Italy from July 2012 to June 2013 and the testing area was Italian southern Apennines, Po Plain, they found that the false positive rate was lesser than 33% while the missing rate is up to 67% (Tramutoli, Corrado et al., 2015).Eleftheriou et al. studied the earthquakes occurred in Greece 5 in period 2004-2013 with use of TIR images acquired by MSG/SEVIRI, more than 93% of all identified TIR anomalies occurred in the prefixed space-time window around time and location of occurrence of earthquakes (with M>4) and the overall false positive rate is <7% (Eleftheriou, Filizzola et al., 2016).It seems that RST is an effect method to extract TIR anomalies before earthquakes, but there is no such study for the mainland of China 10 However, some researchers have proved that some single earthquake results are unreliable.Some so-called TIR anomalies are caused by weather anomalies, which are not related to earthquakes.An instance is that Matthew et al. studied the Gujarat (India) earthquake of 2001, and he found that the previous study, which indicated there was TIR anomalies before the earthquake, was not reliable.They concluded that there was no robust evidence for the existence of LST 15 anomalies prior to the 2001 Gujarat earthquake and cloud covering was one possible cause for the anomalies (Blackett, Wooster et al., 2011).So, a rigorous rule of thermal anomaly judgement and long period statistical analysis are necessary.
In this paper, RST will be applied for the area ( 27 ° − 37 °, 97 ° − 107 ° ) to study a mountainous area of China.The long-term analysis (from Sep.2002 to Mar. 2018) will conform the 20 correspondence between TIR anomalies and earthquakes.Based on the statistical results, the earthquake prediction ability of RST will be evaluated in this paper.

Study area 25
Fig. 1 The distribution of faults in East-Southern Gansu province and its neighbor regions.
East-southern Gansu province and its neighbor regions are selected to study the correlation of the study area is 27 °  37 °, 97 °  107 °.The region is located at the junction of Gansu, Qinghai and Sichuan, it is also the intersection of the North of North-South seismic belt and Kuma seismic belt, moreover, structures in this area are complex and strong earthquakes are frequent (Yang, Zhang et al., 2002).The area is on the eastern edge of Tibet Plateau, belonging to the upper part of the rhombic block in the eastern-southern Gansu province.The Xian River fault, Longmen Shan 5 Fault and Anning River fault are connected here, and the shape of the structure is like a 'Y'.This kind of geomorphology contains abundant plate tectonics, and the Longmen Shan fault in the direction of NNE becomes the steep slope of the Southeast Sichuan Basin and the erosion plateau in the northwest of the study area.10 3. Data and method 3.1 Data introduction MODIS data is used to calculate the TIR anomalies, and the earthquake information got from the China Seismic Information (http://www.csi.ac.cn/) is used for statistical study in this paper.
The MODIS instrument is operating on both Terra and Aqua Spacecrafts.It has a viewing swath 15 width of 2330 km and views the entire surface of the Earth every one to two days.Its detectors measure 36 spectral bands between 0.405 and 14.385 μm, and it acquires data at three spatial resolutions -250m, 500m, 1000m.In this study, nighttime MODIS Land Surface Temperature daily data (MYD11C1) is used for the extraction of TIR anomalies.Because LSTs data are susceptible to solar radiation during the daytime, the nighttime data is selected.The LSTs data is retrieved at 20 5600m by the generalized split-window algorithm.In the day/night algorithm, daytime and nighttime LSTs are retrieved from pairs of day and night MODIS observations in seven TIR bands.Moreover, the daily nighttime Cloud Mask data (MYD35L2) is used for excluding the LST data covered by the cloud.The resolution of Cloud Mask data is 250m and 1000m, so the resolution have to be downscaled to correspond with the LST data.25 The China Seismic Information has recorded earthquakes of whole China since1965.The total number of the recorded earthquakes is 898035.Earthquake caused by block movement and crust compression is a kind of severe geological movement, earthquakes are the instantaneous bursts of accumulated energy, and they may lead to large area of TIR anomalies.The earthquake happens out of the study area will also cause the TIR anomalies near the boundaries in study area, so, for the 30 previous analysis, the earthquakes happened in 25 °  40 °, 95 °  110 ° will be used to study the TIR anomalies in 27 °  37 °, 97 °  107 °.But earthquakes caused by ground subsidence, human factors and so on will not cause TIR anomalies, so the earthquakes with depth=0 should be excluded.Moreover, Tronin indicated that the anomaly was sensitive to crustal earthquakes with a magnitude more than 4.7 and for distance of up to 1000km (Tronin, Hayakawa 35 et al., 2002).So, we select the earthquake with M ≥ 3.5 and depth >0 happened in the area of 25 °  40 °, 95 °  110 ° for study, and after screening, the total number of earthquakes satisfied conditions is 3615.

RST Methodology 40
The RST approach is based on a multi-temporal analysis of historical data set of satellite observations acquired in similar observational conditions (Eleftheriou, Filizzola et al., 2016).Because surface environment is relatively constant, the location of high temperature and low temperature is also relatively steady.Along with the change of time, infrared brightness temperature will change, but it changes in slow speed and small scale with obvious seasonal characteristics.Put aside influence of meteorological factors and earthquake thermal infrared anomaly, the brightness temperature in the same area and same period of the year has strong stability and regularity.So, the basic theory of RST is that the background field is constructed to extract the thermal anomalies, and the mean and variance of the land surface temperature are used to measure t he degree of thermal 5 infrared anomalies.
This method consists of three main steps, it is as follows:  Pre-processing RST is to construct a reference, which is regarded to be at normal state under no influence from other factors, and to measure and extract the anomalies at corresponding time.(, ) is LST data 10 in location r at time t.So, the first step is to eliminate the data affected by clouds, and to remove outliers.
 To eliminate the influence of day to day climatological changes or season time drifts, a pre-processing will be conducted for the daily LST data: ∆V(r,t) is difference between value of LST acquired at time t in location r and its spatial average V(t) computed on the investigated area considering only belonging to the same class, while in this study area all the pixels belong to the land class. To build cloud-mask with use of Cloud Mask data (MYD35L2).In order to be sure that only cloud-free radiances contribute to the computation of reference fields, not only those 20 pixels but also the 24 ones in a 5 × 5 box around it (very often belonging cloud edges) have been excluded by the following computations of reference fields (Eleftheriou, Filizzola et al., 2016).
 To build an outlier-mask.25 Apart from the influence of clouds, there are still many other factors will change the LST, for instance, the extreme weather (blizzard, foehn effect and so on), human's activity and forest fire.These factors will lead to rapid and dramatic changes in large area of LST data.And these anomalies should be excluded from the construction of background filed and the extraction of TIR anomalies.30 The detailed steps for calculating A 3 will be shown in next step.

 Computing Reference Fields (RF)
 The  ∆ (, , ∆) is the mean of location r for time series T and variance  ∆ 2 (, , ) is applied at the time  , by using homogeneous historical records collected under the 35 temporal constraint t ∈ T in the past (t<τ), and the   (, , ∆) is the background field.
A(r, T) = A 1 (, t) * A 2 (, t) * A 3 (, t) (4)  To compute the outlier-mask, and this method is aimed to eliminate the abnormal significant data caused by the non-seismic factors.First step is to find reference fields of minima   (, , ) and maxima   (, , ) to be used at the time τ.
∆  (, , ) ≡ min{∆(,  1 ), ⋯ , ∆(,   )} ∀ ∈ , A(r, T) = 1 (7) ∆  (, , ) ≡ MAX{∆(,  1 ),⋯ , ∆(,   )} ∀ ∈ , A(r, T) = 1 (8) 5Where a modified outlier-mask A 3 (, t) = A(r, T) • (, , ) have been introduced in order to avoid contributions from accidental minima or maxima.And because the blizzard, forest fire and the large area of clouds usually cause the abnormal increase or decrease in LST with a magnitude that far bigger than the change caused by earthquakes.These data should not be included for the calculation of reference field.As shown by ALIANO et al. 10 and GENZANO et al., spatial distribution of clouds, over a thermal heterogeneous scene, can significantly change the value of ∆ in the cloud-free pixels (Saraf, Rawat et al., 2009).The large area of clouds will bring a cold spatial average effect to the computation of reference fields, so when (, ) <   − 2 *   (here,   is the temporal average and the   is its standard), the value of these pixels will be excluded (Eleftheriou, Filizzola 15 et al., 2016).Moreover, even if not producing a cold spatial average effect, an extended cloud coverage can determine values of V(t) and then of the considered signal ∆V(r, t) scarcely representative of the actual conditions of cloud-free pixels, so when the cloudy fraction of land portion of the scene is > 80%, then all the pixels have to be excluded from the computation of reference fields (Eleftheriou, Filizzola et al., 2016).20 The term δ 1 (r,t, τ)  δ 2 (r, t, τ) are defined by the expression: For the outliers caused by the forest fire or blizzards and other factors, δ 3 is used to remove them (with k >=2, in this study the k is set to be 4), and its expression is as follows: 25 The process should be iterated until no further exclusions, are determined by the use of the latest determination of δ (Tramutoli, 1998).This process should be paid more attention, because in the past papers, this process is always ignored.[ ( , ) ( , , ) ( , , ) ( , , ) is affected by cloud.In the results, it can be easily concluded that some area in certain time will lack of data, and for these situations we have to implement a special value to inform that this data is affected by clouds and it should 5 not be analyzed in the following part.

Identification of TIR anomalies
After the calculation of K(r, τ, T), the next step is to determine whether it is the TIR anomalies and to correlate the TIR anomalies with earthquakes.In this paper, the K(r, τ, T) bigger than the 10 threshold means that there is a TIR anomaly, moreover, other conditions will be applied to conform the correlation between the TIR anomalies and earthquakes.For K(r, τ, T) and (, ), only these cases obey the following conditions, can be concluded that the TIR anomaly is related to (, ): In ELEFTHERIOU's study, the threshold was set to 4, however, from the point of view of mathematical statistics, when the value is greater than 15 two times the standard deviation, it has already belonged to the abnormal category, so in this study, the threshold is set to 2.
2) The (, ) is surely not be blocked by clouds or affected by other factors.
3) Spatial persistence The TIR anomalies are gathering together, and they are not isolated being part of a group covering at least 150 2 within an area of 1° * 1° (400 pixels in 20 the images).4) Temporal persistence There are at least one more TIR anomaly appearing after the first TIR anomaly in 7 days.5) The TIR anomalies appears 30 days before or 15 days after the (, ) (Eleftheriou, Filizzola et al., 2016) .25 6) The shortest distance for one point in the TIR anomalies group to the epicenter of (, ) is less than   = 10 0.43 .The cases that the TIR anomalies satisfy conditions 1), 2), 3) and 4) but do not satisfy at least one of 5) and 5) mean that there are TIR anomalies but no corresponding earthquake.And other kinds of cases are of no TIR anomalies.30

Results and analysis
A comprehensive statistical analysis and the TIR extraction results are shown in this chapter.In chapter 4.1, a statistical analysis is conducted to describe to the basic seismological conditions in the study area.And the statistical results for the correlation between earthquakes and TIR anomalies 35 are shown and analyzed in chapter 4.2.Finally, an analysis of earthquake prediction ability for RST will be shown in chapter 4.3.

A basic statistical analysis for earthquakes in study area
Before studying the correlation between TIR anomalies and earthquakes, a simple analysis for 40   And another evidence is shown in the Table .1,we have counted the earthquake with M ≥ 5.0 in period B (from 2002.09 and 2007.12).There are 229 earthquakes with M ≥ 5.0, but the total number in period A is 24, which is accounted for 10.48%, while the duration of period A is accounted for 33.87% of the total time (period A + period B).Moreover, there is no earthquakes with M ≥ 6.5 5 in period A, but there are 14 earthquakes with M ≥ 6.5 in period B. All these evidence tell us that the seismic activity in period B is much more violent and frequent.earthquakes with M≥ 6, the reason is that earthquakes usually occur in faults of active geological movements.
The purpose to study the temporal and spatial characteristic of earthquakes is for a general understanding of the seismic activities in study area.But there is another important reason which is to avoid the large accumulation of earthquakes in a short time and a small area, and that will lead to 10 the same TIR anomaly corresponds to too many earthquakes, this phenomenon will excessively increase the statistical results in the latter part.It is found that there are about 233 earthquakes after the 12 th May 2008 Ms8.0 Wenchuan earthquake, and the locations of these earthquakes are close to epicenter of Wenchuan earthquake.In chapter 4.2, the statistical analysis will be divided into two part, one is with the earthquakes in period C (from 2008.04 to 2008.07), and the other is without the 15 period C.

The statistical analysis for the correlation between TIR anomalies and earthquakes
In this section, the TIR anomalies are extracted and a statistical analysis aiming at studying the correlation between TIR anomalies and earthquakes with M≥ 4 is also been conducted.The way 20 to judge the TIR anomalies is strictly conformed the rules mentioned in chapter 3.3.It can be concluded from the Fig. 5 that 35 TIR anomalies correspond to earthquakes, while the other 26 do not, which are rows 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 14, 15, 16, 17, 19, 20, 21, 22, 23, 35, 41, 43, 5 44,49, 51   However, as it has been mentioned in section 4.1 that period C may increase the total number of TIR anomalies and the correlation rate a lot.So, the experiment without period C is also been conducted, and the number of TIR anomalies is 60, while the correlation rate is 56.7%, both of 15 which are roughly the same as the former result.Theoretically, in period C, the frequency of earthquakes is high and the magnitude of earthquakes is large, which should generate a lot of TIR anomalies and have good correspondences with earthquakes.But there is only one TIR anomaly corresponding to 5 earthquakes, Fig. 7 may indicate the reason.It indicates that earthquakes are gathering in several faults, but the spatial location of these faults are always blocked by the clouds, 20 and the percentage is bigger than 90%.With the long time and large area of clouds blocking, the TIR anomalies caused by the earthquakes in period C may cannot be extracted by RST.
For a more comprehensive analysis for Fig. 5, there are 23 TIR anomalies in period A with 19 corresponding to no earthquakes, while there are 38 TIR anomalies in period B with only seven corresponding to no earthquakes and the correlation rate reaches 82% which is much higher than 25 57%.Fig. 2 and Fig. 8 can explain this phenomenon, in period A, the magnitudes of earthquake intensity are small, the frequency is low, and nearly half of the earthquake occurs in the cloudy region or its adjacent area, so the earthquakes are difficult to correspond to the extracted anomalies, and some potential anomalies may have not been extracted because of clouds.As for the results in period B, the frequency and magnitude of earthquakes in the sparsely clouds areas increase a lot, 30 thermal anomalies are more likely to be extracted and TIR anomalies are more likely to correspond to earthquakes.FN: False Negative, the total number of earthquakes corresponding to no TIR anomalies.PPV: The rate of TIR anomalies which correspond to earthquakes to the total TIR anomalies.FPR: The rate of TIR anomalies which correspond to no earthquakes to the total TIR anomalies.
TPR: The rate of earthquakes which correspond to TIR anomalies to the total earthquakes.FNR: The rate of earthquakes which correspond to no TIR anomalies to the total earthquakes.5 For a more accurate understanding of the 8 parameters, an example is shown in Table 2.The example studies the earthquakes with M ≥ 5.0, and result indicates that there are 61 (TP1+FP) TIR 10 anomalies appearing in study duration, and 19 (TP1) of them correspond to earthquakes while the other 42 (FP) do not, the probability of exact correspondence between the TIR anomalies and earthquakes is 31.1% (PPV) while the probability of no correspondence is 68.9% (FPR).Moreover, there are 250 (TP2+FN) earthquakes with M ≥ 5.0 in the study area, and 32 (TP2) of them correspond to TIR anomalies while the other 218 (FN) correspond to no TIR anomalies, the 15 probability of exact correspondence between the earthquakes and TIR anomalies is 12.8% (TPR) while the probability of no correspondence is 87.2% (FNR).We have calculated the earthquakes with M≥ m (m={3.5, 3.6, 3.7,…,7.8,7.9, 8.0}).And both the experiments with and without period C are calculated, the results show that they are roughly the same, so in this section, we only discuss the result with period C. 20  As it is shown in Fig. 9, PPV declines as the magnitude increases, while FP obviously increases, this phenomenon indicates that with the increase of the magnitude, the number of TIR anomalies corresponding to the earthquakes decreases.TPR and FN decrease steadily, because the total number of earthquake samples is decreasing.5But, what needs more attention is the rates (PPV, TPR).Firstly, a general perceptual analysis shows that PPV decreases steadily with the increasing of M, while TPR is increasing when M≤ 6.6 and M = 7.2, 7.3, the maximum of TPR is 50% when M=6.6 and 7.3, and the TPR decreases when M= 6.7~7.1.When M is 3.5 and 4.0, PPV is 65.6% and 57.4%, it means that when there is a TIR anomaly appearing, there is a possibility of 65.6% (57.4%) that earthquakes with M≥3.5 (4.0) 10 will happen.When the M is 5, 6, 7, the E is 31.1%,19.7%, 3.3%, and these are much lower than PPV of M=3.5 and 4.0.It is concluded from the change of PPV curve that where there is a TIR anomaly, there will be more than 50% possibility of an earthquake with M≥3.5 (4.0) in study area.However, this does not mean that when there is a TIR anomaly, there will be strong earthquakes with M≥5.0 in the study area.On the contrary, the probability of earthquakes with high magnitude 15 is still not high.
The curve of TPR tells the probability of a TIR anomaly when earthquakes occur.When M=3.5 the TPR is 3.1%, and with the increasing of M, TPR increases steadily, but it keeps in a low state when M ∈ [3. 5, 5.4] and the TPR are lower than 20%.The results show that the earthquakes with low magnitude have relatively low possibility (less than 20%) to correspond to TIR anomalies, 20 however, the earthquakes with M ≥ 6.0, which are very destructive, have a relative high correspondence.The high correspondence is of great significance to earthquake prediction.It tells us that a considerable number of the destructive earthquake is more likely to be predicted by this method.
According to two results, it can be concluded that, when a TIR anomaly occurs, there is a 57.3% 25 possibility of an earthquake with M≥4.0, and when there is a strong earthquake with M≥6.0, more than 1/3 of the earthquakes correspond to thermal anomalies.Most TIR anomalies correspond to the earthquakes with M≥4.0, however when M≥6.0, the PPV is relatively low, and that means the false alarm rate for strong earthquakes is high.For TPR, it is increasing with the magnitude, and when M=6.6 and 7.3 it is 50%.It can be concluded from the curve TPR that the greater the 30 magnitude of an earthquake is, the more likely it is to be predicted by this method.But, in fact the PPV and TPR are low, or in other words the FPR and FNR, which are negative for the prediction ability of RST, are really high.All in all, the false alarm rate for M≥4.0 is 42.6%, with the increasing of M the FPR will become higher, and the missing rate for M≥4.0 is 96%, it seems that when M<5.5, there is no obvious correlation between TIR anomalies and earthquakes, though TPR 35 will increase when M is becoming bigger, while the maximum of TPR is still 50%, which is also an unsatisfied value.So, the prediction ability of RST in Sichuan area is limited.

Discussion
To compare the results with the previous similar studies, the summary of the four similar studies 40 is shown in Table .3.It is obvious that PPV of this paper is relatively lower than others, so it is important to verify its actual added value in comparison with a random alarm function (see for instance Eleftheriou).The detailed method is accessible in chapter 3.4 of (Eleftheriou, Filizzola et al., 2016), and the result is shown in Fig. 10.When M≥ 3.5, the point is at the upper part of the random guess, which means that there is no obvious correlation between TIR anomalies and earthquakes with M≥ 3.5, it seems that it is of a casual correlation.When M≥ 4.0 and M≥ 4.5 these two points are still very close to the line, but it is at the lower part, and that means that a noncasual correlation actually exists among the extracted TIR anomalies and earthquakes (M≥ 4.0 and 5 ≥ 4.5).However the correlation is not strong.The result in this paper is different from Eleftheriou's study, in her study, there is a much more obvious and strong correlation between the TIR anomalies and earthquakes.The cause may be that, as it being shown in Fig. 8, east and southeast corner of the study area is always blocked by the clouds, and total time taken up is over 90%, at the same time there are many earthquakes happening in this area, and these earthquakes cannot be related to TIR 10 anomalies because of the lack of data, and they are i nevitably counted into FNR, which is v in Molchan error analysis, making the correlation weaker.For the M≥ 5.0, 5.5 and 6.0, the points are obvious under the random guess, with the increasing of M, the non-causal correlation is becoming stronger.
Tronin indicated that the anomaly was sensitive to crustal earthquakes with a magnitude more 15 than 4.7 and for distance of up to 1000km (Tronin, Hayakawa et al., 2002).But in this study, when M=4.7, the TPR is 9.6% and the FNR is 90.4%, and the point in Fig. 10, when M≥ 4.5, is very close to the random guess, so the statistical result does not support the opinion that TIR anomaly is sensitive to the earthquakes with M ≥ 4.7.When M≥ 5.9, the earthquakes seem to be sensitive to TIR anomalies according to the Table .3.The reason for this failure to conform to previous 20 conclusion may be the regional structure and geological movement, cloud cover, the effectiveness of method extracting TIR anomalies and so on.Further study is needed.We have counted the total number of TIR anomalies and number of FP in each month of two studies.And it is found that the TIR anomalies are gathering in November, September, January and February in Eleftheriou's study, while in this paper, TIR anomalies are gathering in November, September and January.A line which means the percentage of area not blocked by clouds in Sichuan 30 area is also presented in Fig. 11.In this paper, there is a significant positive proportional correlation between the number of TIR anomalies and the area not blocked by clouds.When the percentage is high, the number of TIR anomalies is high, and when the percentage is low, the number of TIR anomalies is low.So there might be many TIR anomalies related to earthquakes, but they are blocked by the clouds and cannot be extracted.However, what is the true cause, blocked by clouds, caused 35 by seasonal weather or else?It is remained to be solved.Moreover, there is another interesting phenomenon that many TIR anomalies corresponding to no earthquakes are gathering in November and September, and both months are not cloudy and are in cold weather.So, why many FP are gathering in these two months is also remained to be studied.anomalies in Eleftheriou's study, 'FP in E' is the number of TIR anomalies corresponding to no earthquakes in Eleftheriou's study.

Conclusion
(1) 18 years statistical analysis for the correlation between earthquakes and TIR anomalies shows 5 that 57.4% TIR anomalies correspond to the earthquakes with M ≥ 4.0 in Sichuan area, and the bigger the M is, the more likely the earthquakes will correspond to TIR anomalies.The low PPV and TPR may be due to the large areas in the study region are covered by clouds all year round.
(2) The low PPV and TPR suggest that the earthquakes prediction ability in Sichuan area with use 10 of RST is limited.For the strong earthquakes with M ≥ 6.0, though the false alarm rate is high, the missing rate is relatively low.RST applied to the study area and has certain prediction ability for strong earthquakes.
(3) There is no obvious correlation between the earthquakes with M < 5.0 and the TIR anomalies extracted by RST in Sichuan area.However, the underlying causes of this situation need to be 15 further studied.(4) RST put up in this paper or in Eleftheriou's study is still seriously affected by the clouds and the seasons.It is necessary to improve and optimize algorithms and statistical methods to exclude the influence of clouds and seasons.20 Algeria, Iran, China, Pakistan and Indonesia through NOAA-AVHRR, Terra/Aqua-MODIS and Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-214Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 August 2018 c Author(s) 2018.CC BY 4.0 License.

35 
Change-detection step  To compute the Index of change of the Environment (ALICE), and the bigger the absolute value is, the more obvious the anomaly is ( , , ) V r T    is the ALICE of location r at Nat. Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-214Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 August 2018 c Author(s) 2018.CC BY 4.0 License.time τ, which is belonged to the time series T.
Nat. Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-214Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 16 August 2018 c Author(s) 2018.CC BY 4.0 License. the temporal and spatial characteristic of earthquakes is conducted.Firstly, the temporal distribution shows that the seismicity in 2008 is the most active from 2002 to 2018, and the seismicity after 2008 is obviously more frequent and violent than it before 2008.The bottom of Fig.2 indicates that there are 3615 earthquakes in the study area, 3.5 ≤ M ≤ 4 is 2262, 4 ≤ M ≤ 5 is 1124, 5 ≤ M ≤ 6 is 198, 6 ≤ M ≤ 7 is 26 and 7 ≤ M ≤ 8 is 5.So the 5 study area is a region with severe earthquake activity.And, as it is shown in upper of Fig.2, the average frequency in period A (from 2002.09 to 2007.12) is about 78.But total number of earthquakes in 2008 increases to 981 including the May 12 2008 Ms8.0 Wenchuan Earthquake, which is the most serious earthquakes in China in recent years.Though, the frequency decreases a lot after 2008 (the average frequency in this period is 243), it is still much higher than the frequency 10 in period A. The temporal distribution tells that the seismic activity before 2008 is relatively weak, but in 2008, the seismic activity is extremely intense and reached the peak value, after 2008 the seismicity in this area still maintains a relatively strong state.

Fig. 2
Fig.2 The temporal distribution from 2002.09-2018.03 of earthquakes with M ≥ 3.5 in study area, 15 and the distribution of seismic frequency along with the magnitude of the earthquakes.Table.1The catalog of earthquakes with  ≥ . before year 2008.

Fig. 3
Fig.3 the spatial distribution of earthquakes in study area, and the brown rectangle refers to the study area (27 °  37 °,97 °  107 °).The reason showing the earthquakes out of study area is that 10 the earthquakes closed to the study area may also cause TIR anomalies in study area.

Fig. 3
Fig.3 shows the spatial distribution of earthquakes.The result indicates that the earthquakes mainly gather in the west and center of study area ( 25 °N to 40 °N, 95 °E to 110 °E) where are mountainous and the earthquakes are mainly aggregated along faults, while the spatial distribution in the east and Sichuan Basin is much more sparse.There is a clustering phenomenon centered on 5

Fig. 4
Fig.4 Two instance for the correlation between TIR anomalies and earthquakes, and the left is the TIR anomaly in 2006.12.29 corresponding to two earthquakes and the right is the TIR anomaly in 2010-10-22 corresponding to three earthquakes.25

Fig. 6
Fig.6 Distribution of TIR anomalies with respect to the earthquakes occurrence for different class 10 of magnitude.

Fig. 7 5 Fig. 8
Fig.7 The distribution of earthquakes in period C and the frequency of each pixel being blocked by cloud in period C, and the higher values mean that the pixels are more frequent blocked by clouds.

Fig. 9
Fig.9 The statistical result of earthquakes including the period C (from 2008.04 to 2008.07), and the curve TP1, FP, TP2, FN, PPV and TPR are corresponding to the example in Table 2. Table.2The detailed results of earthquakes including the period C (from 2008.04 to 2008.07).

5
Fig.11 The monthly average percentage of area not blocked by clouds in study region in Zhang's 10 study, or in other words, in this study.And the bar chart is the numbers of different kind of TIR anomalies.'Total in Z' is the total number of TIR anomalies in this study, 'FP in Z' is the number of TIR anomalies corresponding to no earthquakes in this study.'Total in E' is the total number of TIR t, τ) computed by using an iterative k-clipping technique which starts by computing (, , ) on the base of the first determination of   (, , ) and   2 (, , ), continues by updating their values by using only space/time
and 52.The correlation rate is 57.4%.And as it is shown in Fig 6 that most TIR anomalies appear before earthquakes.

Table . 3
A general summary for research studying the statistical correlation between TIR anomalies and earthquakes.25