A New Algorithm for the Characterization of Thermal Infrared Anomalies in Tectonic Activities

: The monitoring of earthquake events is a very important and challenging task. Remote sensing technology has been found to strengthen the monitoring abilities of the Earth’s surface at a macroscopic scale. Therefore, it has proven to be very helpful in the exploration of some important anomalies, which cannot be seen in a small scope. Previously, thermal infrared (TIR) anomalies have been widely regarded as indications of early warnings for earthquake events. At the present time, some classic algorithms exist, which have been developed to extract TIR anomaly signals before the onset of large earthquakes. In this research study, with the aim of addressing some of the deﬁciencies of the classic algorithm, which is currently used for noise ﬁltering during the process of extracting tectonic TIR anomalies signals, a novel TTIA (tectonic thermal infrared anomalies) algorithm was proposed to characterize earthquake TIR anomalies using the Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperature dataset (MOD11A2). Then, for the purpose of determining the rule of the TIR anomalies prior to large earthquake events, the Qinghai-Tibet Plateau in China was chosen as the study area. It is known that tectonic movements are very active in the study area, and major earthquakes often occur. The following conclusions were obtained from the experimental results of this study: (1) The TIR anomalies extracted using the proposed TTIA method showed a very obvious spatial distribution characteristic along the tectonic faults, which indicated that the proposed algorithm had distinctive advantages in removing or weakening the disturbances of the atectonic TIR anomalies signals; (2) The seismogenic zone was observed to be a more effective observation scale for assisting in the deeper understanding and investigations of the mid- and short-term seismogenic and crust stress change processes; (3) The movement trace of the centroids of the TIR anomalies on the Tibetan Plateau three years prior to earthquake events contributed to improved judgments of dangerous regions where major earthquakes may occur in the future.


Introduction
Earthquake events bring enormous disasters to human society. Therefore, the monitoring of earthquakes is a huge and challenging task with great significance. Traditional earthquake monitoring techniques can only observe the systems of the Earth at very local positions (i.e., deformation observation and underground water temperature and level, outflow, and hydrogen ion exponent (pH), etc. from the ground site) [1][2][3][4], which limits the observations of environments with seismic activities at a macroscopic level. Meanwhile, remote sensing technology has been found to strengthen the monitoring abilities of the Earth's surface at a macroscopic level [5][6][7]. Therefore, these new types of technology have opened up a new era for earthquake monitoring by simultaneously obtaining large amounts of information concerning the dynamic features of the Earth's crust and seismic activities [8][9][10][11][12][13][14].
Enhanced thermal infrared (TIR) emissions from the Earth's surface preceding an earthquake, which are often perceivable by remote sensors, can be referred to as a thermal anomaly [8,9,[15][16][17][18][19][20][21][22][23]. TIR anomalies have been widely reported [16][17][18][19][20] to occur prior to some of the major earthquake events in the past 20 years using satellite data, among which Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) data are the most widely used [3,5,[7][8][9]. Some researchers [2,8,21,43] have pointed out that the TIR anomalies are associated with fault systems or main tectonic regions. These observations have been regarded as important potential earthquake indicators, and have assisted in the predictions of earthquake occurrences. This strategy has attracted a great deal of attention from researchers during recent years, and has contributed to intense discussions on the mechanisms of precursory anomalies [3,[44][45][46][47][48][49], which can be attributed to degassing from rocks under stress, and the participation of ground water have been propounded as a possible cause for generation of TIR anomalies or to p-hole activation in the stressed rock volume and their further recombination at the rock-air interface and frictional heating on fault surfaces as a result of seismicity [9,23,[50][51][52][53][54][55][56]. The recent development of the Lithosphere-Atmosphere-Ionosphere (LAI) coupling model and experimental data of remote sensing satellites on thermal anomalies before major strong earthquakes have described the probable physical basis behind this land surface temperature (LST) anomaly, namely that the radon emanations in the area of earthquake preparation can possibly produce variations of the air temperature and relative humidity [2]. The main physical process responsible for the observed variations maybe is a result of water vapor condensation on ions produced as a result of air ionization by energetic α-particles emitted by 222 Rn [2]. Other intense discussions mainly focuses on explorations regarding spatial-temporal evolution patterns [29,32,57], and the development of algorithms for the identification of TIR anomalies [35,[58][59][60][61][62][63][64][65][66][67][68][69][70][71][72], in which the Robust Satellite data analysis Technique (RST) proposed by Tramutoli in 2007 is a very well-known TIR anomalies algorithm [22] that is widely used to observe the anomalies prior to the major earthquake events [9,21,22,25,31,33,42,43,62]. It has become apparent that the proper presentation of the precursory anomalies is very crucial to the predications of earthquake events. Therefore, the determination of methods for the accurate extraction of the TIR anomalies, which are known to be highly related to tectonic activities, has become imperative. In actual situations, it has been found to be very challenging to extract tectonic TIR anomalies from the entire information of an area, which may include both the tectonic and atectonic heat information. This is due to the fact that the tectonic TIR anomalies caused by the activities of the Earth's crust tends to make up a very small percentage within the entire heat information data. In this study, by reviewing the development of the TIR anomalies algorithms, from the early stage TIR anomalies method characterized by visual interpretation, to the more recent sophisticated algorithms, it was possible to approximately divide the algorithms into three categories as follows: (1) The algorithms based on differential analysis [58][59][60], such as the brightness temperature difference method. In this method, the anomalies are determined by calculating the differences between of the brightness temperature values prior to and after an earthquake event [58]; or calculating the differences between the inner and outer tectonic faults [59,60]. These algorithms are simple. However, they cannot efficiently eliminate the impacts of short-term The Tibetan Plateau is a very important seismically active region which is known for its high incidence level of earthquakes due to its special tectonic region, i.e., locating the convergence region where the Indian Plate extrudes the Eurasian plate from south to north. Most of the tectonic activity occurs on the boundaries of the blocks [75,76]. Earthquakes are the results of abrupt releases of accumulated strain energy that reaches the threshold of strength of the earth's crust. The boundaries of tectonic blocks are the locations of the most discontinuous deformation and highest gradient of stress accumulation, thus they are the most likely places for strain energy accumulation and releases, and in turn, devastating earthquakes, to occur [75,76]. Previous research indicated that the mechanisms of precursory anomalies can be attributed to degassing from rocks under stress, and the participation of ground water or to p-hole activation in the stressed rock volume and their further recombination at the rock-air interface, especially at the boundaries of active tectonic zone [3,9,23,37,[44][45][46][47][48][49][50][51][52][53][54][55][56], thus more TIR anomalies occurred along the tectonic zone [2,8,21,43]. It is wellknown that the Tibetan Plateau is made up of the several important second level active tectonic blocks, detailed in Figure 1. Therefore, the area has been hit by many major earthquakes. Only within this territorial region of China, there were nearly 50 shallow earthquakes (focal depth not more than 20 km) with magnitudes greater than Ms 6.0 during the period examined in this study (2003 to 2015, illustrated in Figure 2). Among these, the earthquakes with magnitudes larger than  The Tibetan Plateau is a very important seismically active region which is known for its high incidence level of earthquakes due to its special tectonic region, i.e., locating the convergence region where the Indian Plate extrudes the Eurasian plate from south to north. Most of the tectonic activity occurs on the boundaries of the blocks [75,76]. Earthquakes are the results of abrupt releases of accumulated strain energy that reaches the threshold of strength of the earth's crust. The boundaries of tectonic blocks are the locations of the most discontinuous deformation and highest gradient of stress accumulation, thus they are the most likely places for strain energy accumulation and releases, and in turn, devastating earthquakes, to occur [75,76]. Previous research indicated that the mechanisms of precursory anomalies can be attributed to degassing from rocks under stress, and the participation of ground water or to p-hole activation in the stressed rock volume and their further recombination at the rock-air interface, especially at the boundaries of active tectonic zone [3,9,23,37,[44][45][46][47][48][49][50][51][52][53][54][55][56], thus more TIR anomalies occurred along the tectonic zone [2,8,21,43]. It is well-known that the Tibetan Plateau is made up of the several important second level active tectonic blocks, detailed in Figure 1. Therefore, the area has been hit by many major earthquakes. Only within this territorial region of China, there were nearly 50 shallow earthquakes (focal depth not more than 20 km) with magnitudes greater than Ms 6.0 during the period examined in this study (2003 to 2015, illustrated in Figure 2). Among these, the earthquakes with magnitudes larger than   Figure 1) [77]. The main earthquake event was followed by thousands of aftershocks. Among these, 53 reached or exceeded Mw 5.0. This major earthquake caused extreme damages that extended almost 300 km along the Longmen Shan and its adjacent area. A large number of landslides, debris flows, surface fractures, and dammed lakes were observed. Overall, more than 80,000 fatalities were reported. This catastrophic event was one of the deadliest earthquakes in China during the past few centuries.
On 12 February 2014, at 17:19 local time (09:19 UTC), a significant Ms 7.3 Yutian (Xinjiang, China) earthquake occurred in the Altyn-Haiyuan fault zone of the northwest edge of China's Tibetan Plateau (represented by the yellow square box labelled 1 in Figure 1). The epicenter was located at 36.1°N, 82.5°E, and the depth of the hypocenter of 12 km. This major earthquake event was followed by 739 aftershocks, one of which reached Ms 5.0, and ten reached or exceeded Ms 4.0. As a result of this major disaster, more than 455,000 people became homeless; 68,340 houses were toppled; 11,515 domestic animals were killed; 497 bridges were badly damaged; and 113 landslide events occurred.

Remotely Sensing Data
The LST is an important index for evaluating the heat balance, which represents a collective outcome of the tectonic activities of the deep Earth's crust. Also, the terrain and land cover and meteorological activities are known to be important indexes. Therefore, this study adopted the LST to measure tectonic activities. The data used in this study were the synthetic LST datasets with the average surface temperatures of eight days (MOD11A2 Version 5 (V5) of MODIS) provided by NASA, with a 1 km spatial resolution. This V5 product was a significant improvement over past versions. In this experiment, in order to weaken the influences of solar radiation and local variations due to cloud cover or shadows during the day-time hours, the data observed at night were used. In addition, since the MOD11A2 was one of the mature LST products, there was no need to implement atmosphere corrections [66]. The detailed information regarding the data set is summarized in Table  1. The original images of the MOD11A2 are shown in Figure 3 according to the time sequence. Figure  3 shows that the original LST of the MOD11A2 were noticeably dominated by solar radiation and terrain, and exhibited very distinct seasonal characteristics. Therefore, for extracting out the TIR anomalies closely related to tectonic activities, the influences of atectonic factors must be removed from the original MODIS LST data and weakened.  Figure 1) [77]. The main earthquake event was followed by thousands of aftershocks. Among these, 53 reached or exceeded Mw 5.0. This major earthquake caused extreme damages that extended almost 300 km along the Longmen Shan and its adjacent area. A large number of landslides, debris flows, surface fractures, and dammed lakes were observed. Overall, more than 80,000 fatalities were reported. This catastrophic event was one of the deadliest earthquakes in China during the past few centuries.
On 12 February 2014, at 17:19 local time (09:19 UTC), a significant Ms 7.3 Yutian (Xinjiang, China) earthquake occurred in the Altyn-Haiyuan fault zone of the northwest edge of China's Tibetan Plateau (represented by the yellow square box labelled 1 in Figure 1). The epicenter was located at 36.1 • N, 82.5 • E, and the depth of the hypocenter of 12 km. This major earthquake event was followed by 739 aftershocks, one of which reached Ms 5.0, and ten reached or exceeded Ms 4.0. As a result of this major disaster, more than 455,000 people became homeless; 68,340 houses were toppled; 11,515 domestic animals were killed; 497 bridges were badly damaged; and 113 landslide events occurred.

Remotely Sensing Data
The LST is an important index for evaluating the heat balance, which represents a collective outcome of the tectonic activities of the deep Earth's crust. Also, the terrain and land cover and meteorological activities are known to be important indexes. Therefore, this study adopted the LST to measure tectonic activities. The data used in this study were the synthetic LST datasets with the average surface temperatures of eight days (MOD11A2 Version 5 (V5) of MODIS) provided by NASA, with a 1 km spatial resolution. This V5 product was a significant improvement over past versions. In this experiment, in order to weaken the influences of solar radiation and local variations due to cloud cover or shadows during the day-time hours, the data observed at night were used. In addition, since the MOD11A2 was one of the mature LST products, there was no need to implement atmosphere corrections [66]. The detailed information regarding the data set is summarized in Table 1. The original images of the MOD11A2 are shown in Figure 3 according to the time sequence. Figure 3 shows that the original LST of the MOD11A2 were noticeably dominated by solar radiation and terrain, and exhibited very distinct seasonal characteristics. Therefore, for extracting out the TIR anomalies closely related to tectonic activities, the influences of atectonic factors must be removed from the original MODIS LST data and weakened.

Introduction of the New Algorithm
The objective of the developed algorithm was to extract the tectonic activities information by filtering or weakening the impact of the atectonic factors (solar radiation, atmospheric disturbances, and human activities. It is widely known that solar radiation is the first and foremost factor which affects the land surface temperature field. It has been regarded as a type of stable annual periodic signal that has observable seasonal fluctuating properties. In addition, the influences of air masses and human activities (urban heat islands) should not be ignored. In accordance with the spatiotemporal scale characteristics of the various interference factors, this study proposed a new algorithm (TTIA) to characterize the tectonic TIR anomalies. The process was introduced in detail as follows: Step (1): Constructing the LST background field and calculating the residuals: The annual trend values were extracted from the sequential LST data, such as the temperature background field, using a harmonic analysis fitted curve method. The detailed introduction on this step is shown in Section 2.3.1. The fitting residual error image (shown in Step ③ of Figure 4) filtered the solar radiation influence and climate change. We then applied a 1-D wavelet transform, and by means of deleting the first order high frequency of wavelet transform, the effect of the short-term meteorological factors could be removed. However, it still consisted of disturbances in the atmosphere and human activities.
Step (2): Spatial filtering to weaken the impacts of the atmosphere and human activities (urban heat island): It was considered in this study that the air masses and urban heat islands did indeed have evident influences on the surface temperature field. Therefore, these disturbance factors needed to be removed or weakened. The spatial scales dominated by these two factors were profoundly different. Generally speaking, an air mass can cover hundreds to thousands of km. However, the influence range of the heat islands only measure tens of km [78]. Therefore, a 2-D wavelet transform technique was used in this study to filter the disturbance factors, and extract the tectonic thermal information. A detailed introduction to this is shown in Section 2.3.2. Step

Introduction of the New Algorithm
The objective of the developed algorithm was to extract the tectonic activities information by filtering or weakening the impact of the atectonic factors (solar radiation, atmospheric disturbances, and human activities. It is widely known that solar radiation is the first and foremost factor which affects the land surface temperature field. It has been regarded as a type of stable annual periodic signal that has observable seasonal fluctuating properties. In addition, the influences of air masses and Remote Sens. 2018, 10, 1941 8 of 33 human activities (urban heat islands) should not be ignored. In accordance with the spatio-temporal scale characteristics of the various interference factors, this study proposed a new algorithm (TTIA) to characterize the tectonic TIR anomalies. The process was introduced in detail as follows: Step (1) Constructing the LST background field and calculating the residuals: The annual trend values were extracted from the sequential LST data, such as the temperature background field, using a harmonic analysis fitted curve method. The detailed introduction on this step is shown in Section 2.3.1. The fitting residual error image (shown in Step 3 of Figure 4) filtered the solar radiation influence and climate change. We then applied a 1-D wavelet transform, and by means of deleting the first order high frequency of wavelet transform, the effect of the short-term meteorological factors could be removed. However, it still consisted of disturbances in the atmosphere and human activities.
Step (2) Spatial filtering to weaken the impacts of the atmosphere and human activities (urban heat island): It was considered in this study that the air masses and urban heat islands did indeed have evident influences on the surface temperature field. Therefore, these disturbance factors needed to be removed or weakened. The spatial scales dominated by these two factors were profoundly different. Generally speaking, an air mass can cover hundreds to thousands of km. However, the influence range of the heat islands only measure tens of km [78]. Therefore, a 2-D wavelet transform technique was used in this study to filter the disturbance factors, and extract the tectonic thermal information. A detailed introduction to this is shown in Section 2.3.2.
Step (3) Presenting the tectonic thermal anomalies information by calculating the value image (shown in 5 of Figure 4): Further details of this calculation process are shown in Section 2.3.3.
The above steps were the core parts of this study's proposed new algorithm, and the process is vividly illustrated in Figure 5. The above steps were the core parts of this study's proposed new algorithm, and the process is vividly illustrated in Figure 5.

Construction of the LST Background Field and Calculation of the Residuals
First of all, this study ranked the pixels according to their temporal sequence. Then, a yearly trend value was extracted as the background field by approximating the curve of the LST over a period of 13 years (2003 to 2015). A key point of this process was that all the observed values were regarded as a whole, and a new approximation function was obtained [79][80][81]. A least squares criterion was used to build the new function as follows:  Least squares approximation By assuming that the known function ( ) f x was defined in the space [ , ] a b , then Where    x is a weight function that is not identical to zero. Also, by assuming that the norm of ( ) f x exists, and that the inner product space is   2 , L a b , then the norm of ( ) f x can be obtained using Equation (1) [82] as follows: The nature of the optimal approximation in 2 L is that, for any     should be looked for on the basis space   which is subject to the following: Equation (3) indicated that the square of the error function reached a minimum value after solving the integral. Therefore,    x can be referred to as the optimal square approximation of ( ) f x .

Construction of the LST Background Field and Calculation of the Residuals
First of all, this study ranked the pixels according to their temporal sequence. Then, a yearly trend value was extracted as the background field by approximating the curve of the LST over a period of 13 years (2003 to 2015). A key point of this process was that all the observed values were regarded as a whole, and a new approximation function was obtained [79][80][81]. A least squares criterion was used to build the new function as follows:

Least squares approximation
By assuming that the known function f (x) was defined in the space [a, b], then ρ(x) ≥ 0. Where ρ(x) is a weight function that is not identical to zero. Also, by assuming that the norm of f (x) exists, and that the inner product space is L 2 [a, b], then the norm of f (x) can be obtained using Equation (1) [82] as follows: The nature of the optimal approximation in L 2 is that, for any f (x) ∈ (a, b), function ϕ(x) should be looked for on the basis space M = {ϕ 1 , ϕ 2 , . . . , ϕ n } as follows: which is subject to the following: Equation (3) indicated that the square of the error function reached a minimum value after solving the integral. Therefore, ϕ(x) can be referred to as the optimal square approximation of f (x).

Optimal Fourier approximation
When the function to be approximated is periodic, or the primitive function is too complex, then the primitive function has to be approximated using discretization processing. In this study's case, this step was performed using a Fourier approximation. By letting M be made up of trigonometric functions (sine or cosine functions) [83,84], then ϕ(x) will be the optimal Fourier approximation of ϕ(x).

Extraction of the LST's yearly trend
In this study, by taking a pixel as an example, and building a Fourier approximation function as described in Equation (4), a fitting process for the discrete sequence of the LST values in the pixel could then be performed: where i = 1, 2, . . . , m, and 2m is the number of harmonic; y represents the approximate values; t is the specific moment in a time sequence; a i denotes the amplitude of the i th sine harmonic; b i is the phase of the i th harmonic; c i represents the amplitude of the i th cosine harmonic; ω i is the phase of the i th harmonic; and ε is the fitting error. By considering the annual and diurnal changes of the LST, two types of harmonics with different frequencies were chosen to approximate the time series. One type with a 0 frequency indicated the average of LST, and that there was no phase. Meanwhile, the other type had a period of eight days, which represented the short-term variations of the LST.
As indicated above, the principle of the least square method was considered in the calculation of the variables during the approximation process, as described in Equation (5). The values of a i , b i , and c i could be obtained by computing derivatives as follows:

Tectonic Thermal Signal Extraction by Spatial Two-Dimensional Wavelet Filtering for the Purpose of Weakening the Disturbances of the Atmospheric and Human Activities
Although the LST fitting residuals had filtered out the main disturbances of the solar radiation, it still contained the disturbances resulting from the atmosphere and human activities. It was believed that separating the two types of thermal information from the residual images would allow for improvements in the extraction of the thermal information related to the tectonic activities. After considering the spatial affected scope of the air masses and human activities, this study applied the different features between them to achieve the signal separation using a 2-D wavelet transform.

2-D continuous wavelet transform
The 2-D continuous wavelet transform was defined as follows [7]: Its inverse was given by Equation (8): where a is a scaling transformation parameter; and b is a displacement transformation parameter. When the value varied, the wavelet function changed, which corresponded to the elongation and contraction. That is to say, when a < 1, Ψ(x, y) was elongated; and when a > 1, Ψ(x, y) had contracted. When the value of b changed, the wavelet function exhibited a displacement.

2-D discrete wavelet transform
Due to the computational complexity of the continuous wavelet transform in terms of the continuous values a, a 2-D discrete wavelet transform is often used to conduct the wavelet transformation of a digital image. In this study, the 2-D discrete wavelet transform was defined as follows [83]: where m, n are the integer,a 0 = 1 and b 0 = 0. The corresponding discrete wavelet transform can then be defined as follows: The reconstruction process was the inverse process of the 2-D discrete wavelet image decomposition. The digital images, which had met certain requirements, could be reconstituted through controlling the high and low frequency components of the decomposition process of the wavelet transform.

Tectonic thermal infrared signal extraction based on a 2-D wavelet transform
Except the solar radiation, among the influencing factors that dominate the land surface temperature, atmospheric air mass and human activities are two important and ought not to be missed factors [84]. Generally speaking, a distributed air mass extends up to a large scope at a horizontal direction with a homogeneous temperature and humidity characterization [85]. In terms of the thermal property, an air mass can be divided into a cold air mass and warm air mass [86]. The arrival of an air mass with certain thermodynamic properties can vary the land surface temperature. Therefore, the influence of an air mass on the land surface temperature must be removed [86]. According to the field of meteorology and physical geography, an air mass can cover a range of several hundred to several thousands of km. In this study, in order to facilitate the calculation process, the temperature components with ranges of more than 1000 km were taken as the contribution of the atmospheric air mass [84][85][86][87]. In regard to the disturbances caused by human activities, such as the urban heat island effect, the influence ranges were roughly tens of km (the length of the diameter of a city) [88,89]. Therefore, approximately 30 km was chosen in this study as the influence range of the human activities [84,89]. Then, using a wavelet transform technique, these two types of disturbances could be weakened to a certain degree.
To be more specific, since the spatial resolution of the MODIS LST dataset used in this study was 1 km, and considering the fact that the sampling of the wavelet transform abided by the geometry power of two, the thermal infrared signal was able to extract more information related to the tectonic activities by using the 2-D wavelet transform from the residual images. This was accomplished by weakening the influences of the solar radiation, atmospheric air masses, and human activities. As a result, it could be defined as T tec , and the calculation formula is given in Equation (11): where T L4 represents the LST, which was extracted from the low frequency signal of the 2-D wavelet transform in the 4th scale (corresponding to a 32 km resolution) of the residual images; T L9 is the LST extracted from the low frequency signal of the 2-D wavelet transform in the 9th scale (corresponding to a 1024 km resolution) of the residual images.

"kσ"rule and offset index K tec
In this study, the "kσ"rule was used to determine the signals of the thermal anomalies. The "kσ" principle was proposed initially by Wright in 1884 [90][91][92][93]. This simple and basic principle is as follows: where x and σ are the mean and standard deviation of the result x, respectively. Here the "kσ" rule is introduced to characterize the variation of land surface temperature observations around their mean, where k refers to the zoom multiples of standard deviation σ. When the result x meets the above formula, then it can be defined as an anomaly signal, as illustrated in Figure 6. where represents the LST, which was extracted from the low frequency signal of the 2-D wavelet transform in the 4th scale (corresponding to a 32 km resolution) of the residual images; is the LST extracted from the low frequency signal of the 2-D wavelet transform in the 9th scale (corresponding to a 1024 km resolution) of the residual images.

Expressions of the Tectonic Thermal Infrared Anomalies
 " "rule and offset index In this study, the "  k "rule was used to determine the signals of the thermal anomalies. The "  k "principle was proposed initially by Wright in 1884 [90][91][92][93]. This simple and basic principle is as follows: where x and σ are the mean and standard deviation of the result x, respectively. Here the "  k "rule is introduced to characterize the variation of land surface temperature observations around their mean, where k refers to the zoom multiples of standard deviation σ. When the result x meets the above formula, then it can be defined as an anomaly signal, as illustrated in Figure 6. The effective thermal infrared signal values were achieved using Equation (11). Then, the thermal anomalies were extracted through calculating the offset index according to Equation (13). In fact, is just the embodiment of the k value in the formula (12). The introduction of aims to characterize the variation of land surface temperature observations around their mean. Meanwhile, can also be perceived as zoom multiples of standard deviation ( ).
where denotes the location of pixel on the remote sensing image, with the abscissa and ordinate ( , ); represents the time of the acquisition of the satellite image, with ∈ , where defines the homogeneous domain of the satellite imagery collected in the same time-slot (hour) of the day and period (month) of the year; is the effective thermal infrared signal; and and are the mean value and standard deviation of the in the same time-slot, respectively. It is worth noting that although equation 13 has a similar formation to the RST method, which was proposed by Tramutoli 1998 (RAT) [25], Tramutoli 2005 (RST) [38] and applied by Lisi et al. 2015 [42], who have used in the same way the reference fields of mean and standard deviation; nevertheless, in equation 13, has a distinct meaning from T in the RST method. Here has filtered the influence of air mass and human activities by a 2-D wavelet transformation, therefore , and in the Equation (13) has distinct basis, compared with the and in the RST method.  k Figure 6. Sketch map of the "kσ" rule for detecting anomalies.
The effective thermal infrared signal values T tec were achieved using Equation (11). Then, the thermal anomalies were extracted through calculating the offset index K tec according to Equation (13). In fact, K tec is just the embodiment of the k value in the formula (12). The introduction of K tec aims to characterize the variation of land surface temperature observations around their mean. Meanwhile, K tec can also be perceived as zoom multiples of standard deviation σ T tec (r i ).
where r i denotes the location of pixel i on the remote sensing image, with the abscissa and ordinate (x i , y i ); t represents the time of the acquisition of the satellite image, with t ∈ τ, where τ defines the homogeneous domain of the satellite imagery collected in the same time-slot (hour) of the day and period (month) of the year; T tec is the effective thermal infrared signal; and µ T tec and σ T tec are the mean value and standard deviation of the T tec in the same time-slot, respectively. It is worth noting that although equation 13 has a similar formation to the RST method, which was proposed by Tramutoli 1998 (RAT) [25], Tramutoli 2005 (RST) [38] and applied by Lisi et al. 2015 [42], who have used in the same way the reference fields of mean and standard deviation; nevertheless, in equation 13, T tec has a distinct meaning from T in the RST method. Here T tec has filtered the influence of air mass and human activities by a 2-D wavelet transformation, therefore µ T tec , and σ T tec in the Equation (13) has distinct T tec basis, compared with the µ and σ in the RST method.

Results
In the study results, two important observation scales were adopted to describe the TIR anomalies in tectonic activities based on the proposed method. The first method was the whole Qinghai-Tibet Plateau scale, while the other was the earthquake generating fault zone scale. The Tibet Plateau contributes to observing the TIR in tectonic activities on a large scale, and to understanding the TIR anomalies evolution of the pregnant process of earthquakes before a large earthquake event. It is well known that earthquakes are the manifestations of abrupt huge energy releases that are directly caused by tectonic activities. Since earthquakes are known to be highly related to tectonic activities, it can be deduced that the relevant precursor signals of earthquakes should be evidently strongly correlated with tectonic structures. We chose two large earthquake events that had occurred on the Tibet Plateau, based on the consideration that before the large earthquakes events occurred, the TIR in the tectonic activities often varied significantly and observably, which makes it more convenient for us to observe the spatio-temporal evolution characterization of TIR and validate the new algorithm.
The period in the red typeface in Figure 7 indicates the co-seismic period. The thermal anomaly images shown in Figure 7(1)-(9) were calculated according to the following steps: Step (1) The calculation of the K tec values was performed using Formula (12), which indicated the tectonic thermal anomaly; Step (2) The abnormal periods were identified, which needed to meet the request that K mean ≥ µ + σ. In this study, K mean is the spatial mean of the K tec value image in a period; µ is the mean of the K tec value images in the multi-periods; and σ is the variance of the K tec value images in the multi-periods; Step (3) All of the selected images were overlaid by Step (2), and the average value of the image was calculated; Step (4) The differences in the value image were calculated between the average value image calculated by Step (3), and the K tec mean value image of the aseismic periods (without the occurrence of large earthquakes). Figure 7 illustrates the fact that the spatial form of the TTIA-based TIR anomalies corresponded well in this study to the tectonic zone distributions, which potentially provide evidence that judging the extracted TIR anomalies signals based on the new method is reasonable. This is particularly evident in Figure 7(3). The noticeable characteristics of the spatial distribution of most of the TIR anomalies from east to west, along with the major faults zones of the entire Tibetan Plateau, leads us to believe that the proposed TTIA algorithm was able to reflect tectonic activities well, as it can efficiently filter out the majority of the other atectonic disturbance information originating from the impacts of solar radiation, land cover, atmospheric and human activities. Figure 7 displays the evolution history of the TIR anomalies before and after the earthquake events, which was helpful in understanding the seismogenic stages of major earthquakes. It is well known that every earthquake has a long seismogenic stage, which may range from several decades to hundreds or even thousands of years. Generally speaking, the larger magnitude earthquakes have longer seismogenic stages. In this study, the TIR anomalies evolutions were observed and analyzed, with special emphasis on the period from the three-year pre-seismic stage to the breaking out stage of the earthquakes, in order to grasp the evolutionary characteristics of the TIR anomalies signals during the latter seismogenic stage. Figure 7(1)-(3) illustrate the developing processes of the TIR anomalies from May of 2005 to February of 2008 (prior to the Wenchuan earthquake) over the entire Tibetan Plateau region. The value of the TIR anomalies can be obviously seen as gradually rising during the period from three years to three months prior to the Wenchuan event. Also, in that period the K tec value of the east was clearly higher than that of the west on the plateau. According to the theory of the relationship between the heat and stress in laboratory testing (temperature-rise corresponding to extrusion stress), it was deduced that the surrounding rock in the active tectonic block was increasingly extruded along the direction from east to west. This was confirmed by the observed GPS results. Until February 9th, 2008 (approximately three months before the Wenchuan earthquake), the temperature in the area had been reaching its highest value, as shown in Figure 7(3)). This was the first observed important temperature-rise stage. Then, after passing the stage of the highest value of the TIR anomalies, the temperature was observed to slightly decline, as shown in Figure 7(4). Soon after, the devastating earthquake event occurred, after which the TIR anomalies fell to a relatively low level, as detailed in Figure 7 Figure 9 illustrates the moving tendency of the intensity centroid in the K tec evolution process prior to the Wenchuan earthquake event. The green, yellow, and red points denote 2 to 3 years, 1 to 2 years, and 3 months to 1 year prior to the earthquake, respectively. The centroid apparently had gradually moved from the left to the right of the TIR anomalies image, and toward the epicenter of the Wenchuan earthquake. In reality, this was a significant signal, and important information was obtained that revealed that the eastern region over the Tibetan Plateau seemed to be facing a higher threat than the western region, which will potentially contribute to future estimations of dangerous regions where earthquakes might occur.
Remote Sens. 2017, 9, x FOR PEER REVIEW 16 of 34 Figure 9 illustrates the moving tendency of the intensity centroid in the evolution process prior to the Wenchuan earthquake event. The green, yellow, and red points denote 2 to 3 years, 1 to 2 years, and 3 months to 1 year prior to the earthquake, respectively. The centroid apparently had gradually moved from the left to the right of the TIR anomalies image, and toward the epicenter of the Wenchuan earthquake. In reality, this was a significant signal, and important information was obtained that revealed that the eastern region over the Tibetan Plateau seemed to be facing a higher threat than the western region, which will potentially contribute to future estimations of dangerous regions where earthquakes might occur. In regard to the Yutian earthquake, the evolution trend of the can be seen in Figure 10. As can be seen in the figure, the value of was not only consistently distributed along the trend of the tectonic faults, but also the of the surrounding rock in the tectonic block was observed to steadily rise until the highest value was reached. At that time, the value of the of the western region was higher than that of the eastern region.
Throughout the entire process of the evolution over the Tibetan Plateau in both of the major earthquake cases examined in this study, it was found that prior to the earthquake events, the spatial distributions of were almost entirely arranged along the significant fault zones. Also, the majority were extended along an east-west direction, which confirmed that the proposed TTIA algorithm could accurately reflect the tectonic activities. Therefore, it was concluded that the TTIAbased could be regarded as a precursor indicator for outlining potentially dangerous regions where major earthquakes may occur in the future. In regard to the Yutian earthquake, the evolution trend of the K tec can be seen in Figure 10. As can be seen in the figure, the value of K tec was not only consistently distributed along the trend of the tectonic faults, but also the K tec of the surrounding rock in the tectonic block was observed to steadily rise until the highest value was reached. At that time, the value of the K tec of the western region was higher than that of the eastern region.
Throughout the entire process of the K tec evolution over the Tibetan Plateau in both of the major earthquake cases examined in this study, it was found that prior to the earthquake events, the spatial distributions of K tec were almost entirely arranged along the significant fault zones. Also, the majority were extended along an east-west direction, which confirmed that the proposed TTIA algorithm could accurately reflect the tectonic activities. Therefore, it was concluded that the TTIA-based K tec could be regarded as a precursor indicator for outlining potentially dangerous regions where major earthquakes may occur in the future.
Remote Sens. 2017, 9, x FOR PEER REVIEW 16 of 34 Figure 9 illustrates the moving tendency of the intensity centroid in the evolution process prior to the Wenchuan earthquake event. The green, yellow, and red points denote 2 to 3 years, 1 to 2 years, and 3 months to 1 year prior to the earthquake, respectively. The centroid apparently had gradually moved from the left to the right of the TIR anomalies image, and toward the epicenter of the Wenchuan earthquake. In reality, this was a significant signal, and important information was obtained that revealed that the eastern region over the Tibetan Plateau seemed to be facing a higher threat than the western region, which will potentially contribute to future estimations of dangerous regions where earthquakes might occur. In regard to the Yutian earthquake, the evolution trend of the can be seen in Figure 10. As can be seen in the figure, the value of was not only consistently distributed along the trend of the tectonic faults, but also the of the surrounding rock in the tectonic block was observed to steadily rise until the highest value was reached. At that time, the value of the of the western region was higher than that of the eastern region.
Throughout the entire process of the evolution over the Tibetan Plateau in both of the major earthquake cases examined in this study, it was found that prior to the earthquake events, the spatial distributions of were almost entirely arranged along the significant fault zones. Also, the majority were extended along an east-west direction, which confirmed that the proposed TTIA algorithm could accurately reflect the tectonic activities. Therefore, it was concluded that the TTIAbased could be regarded as a precursor indicator for outlining potentially dangerous regions where major earthquakes may occur in the future. . Figure 10. TIR anomalies spatio-temporal presentation of the Ms 7.3 Yutian earthquake based on the proposed TTIA method on the Tibetan Plateau. Figure 11 shows the TIR anomalies based on the TTIA method in the Altyn-Haiyuan fault zone where the Ms 7.3 Yutian earthquake occurred. It can also be clearly seen that the TIR anomalies are in high correspondence with the tectonic line. Figure 11. TIR anomalies prior to the Yutian earthquake event in different periods. Figure 12 illustrates the moving tendency of the intensity centroid in the evolution process prior to the Yutian earthquake event. The green, yellow, and red points represent the same as in Figure 9. Overall, it was observed that the intensity centroid of the hovered around the western region of the plateau. In particular, in the entire energy of the anomalies with ≥ 1.0 detailed in Figure 12a, the centroid (indicated by the yellow point) almost coincided with the epicenter of the Yutian earthquake. The calculated intensity centroid potentially indicated that the western region of the plateau seemed more dangerous than the eastern region. These findings were very valuable signals, which provided supplementary information for outlining dangerous regions, and served as early warning signs of the earthquake events.  . Figure 10. TIR anomalies spatio-temporal presentation of the Ms 7.3 Yutian earthquake based on the proposed TTIA method on the Tibetan Plateau. Figure 11 shows the TIR anomalies based on the TTIA method in the Altyn-Haiyuan fault zone where the Ms 7.3 Yutian earthquake occurred. It can also be clearly seen that the TIR anomalies are in high correspondence with the tectonic line. Figure 11. TIR anomalies prior to the Yutian earthquake event in different periods. Figure 12 illustrates the moving tendency of the intensity centroid in the evolution process prior to the Yutian earthquake event. The green, yellow, and red points represent the same as in Figure 9. Overall, it was observed that the intensity centroid of the hovered around the western region of the plateau. In particular, in the entire energy of the anomalies with ≥ 1.0 detailed in Figure 12a, the centroid (indicated by the yellow point) almost coincided with the epicenter of the Yutian earthquake. The calculated intensity centroid potentially indicated that the western region of the plateau seemed more dangerous than the eastern region. These findings were very valuable signals, which provided supplementary information for outlining dangerous regions, and served as early warning signs of the earthquake events.  Figure 12 illustrates the moving tendency of the intensity centroid in the K tec evolution process prior to the Yutian earthquake event. The green, yellow, and red points represent the same as in Figure 9. Overall, it was observed that the intensity centroid of the K tec hovered around the western region of the plateau. In particular, in the entire energy of the anomalies with K tec ≥ 1.0 detailed in Figure 12a, the centroid (indicated by the yellow point) almost coincided with the epicenter of the Yutian earthquake. The calculated intensity centroid potentially indicated that the western region of the plateau seemed more dangerous than the eastern region. These findings were very valuable signals, which provided supplementary information for outlining dangerous regions, and served as early warning signs of the earthquake events. In this study result, another observable scale is the earthquake generating fault zone. The spatiotemporal evolution of the TIR anomalies in tectonic activities with the TTIA algorithm in the fault zone scale (Longmenshan fault zone and Altyn fault zone) are shown in Figures 13-16. The following observations can be made: (1) Before the great earthquake, the fluctuations of the value are significantly more volatile than those in an aseismic period, which indicates that the TIR anomalies in tectonic activities in the event year are more active than those in non-event years. (2) In an unperturbed period (when no large earthquake occurs), the value presents low value features in both the Longmenshan fault zone (shown in Figure 14) and Altyn fault zone (shown in Figure 16). Despite the fact that there are also some high or low pixels, their distribution still did not present distinct regularity.
Further observing Figure 13, we can observe the fact that during the period of three to six months before the Wenchuan earthquake event, the surrounding rock of the Longmenshan Fault presented an increasing trend in the , approximately 100 days before the earthquake, and the of the surrounding rock reached its highest value. Then the began to decline. On the contrary, during this period, the fault region went through a cold status before the earthquake. Then around three months before the event, the situation changed. The fault region was going through a hot thermal anomaly state, breaking away from the previous cold anomaly state; i.e., the temperature of the fault zone underwent a change process from cold to hot status. After the earthquake, the value of the in the fault region declined, exhibiting an increasingly discrete state in terms of spatial distribution. A similar phenomenon can be also observed in the Yutian earthquake event which occurred on the Altyn fault on 12 February 2014.
The above mentioned observed phenomenon appears to not be coincidental. In fact, we have investigated the characterization of TIR anomalies in tectonic activities before an earthquake event year via the TTIA method using four earthquake cases of magnitudes greater than Ms 7.0 over the Tibet Plateau since 2003, similar phenomena were observed in each; namely, before the large earthquake, the temperature of the surrounding rock underwent a heating to cooling process, while the temperature of the fault region underwent an inverse process. However, when we investigated earthquakes of magnitude Ms 6.0 or less, we did not observe this kind of very obvious regularity.
The possible explanation for the TIR anomalies in tectonic activities occurring before a large earthquake event is that it is the integration result of the green house air degassing from the Earth's crust and the stress of rocks. On one hand, green house air is able to warm up the air, thereby leading to three phase changes of water vapor, which can produce the variation of the air temperature. On the other hand, the compression and tension of rocks can result in LST variation through a series of complicated geophysical and geochemical processes. The following observations can be made: (1) Before the great earthquake, the fluctuations of the K tec value are significantly more volatile than those in an aseismic period, which indicates that the TIR anomalies in tectonic activities in the event year are more active than those in non-event years. (2) In an unperturbed period (when no large earthquake occurs), the K tec value presents low value features in both the Longmenshan fault zone (shown in Figure 14) and Altyn fault zone (shown in Figure 16). Despite the fact that there are also some high or low K tec pixels, their distribution still did not present distinct regularity.
Further observing Figure 13, we can observe the fact that during the period of three to six months before the Wenchuan earthquake event, the surrounding rock of the Longmenshan Fault presented an increasing trend in the K tec , approximately 100 days before the earthquake, and the K tec of the surrounding rock reached its highest value. Then the K tec began to decline. On the contrary, during this period, the fault region went through a cold status before the earthquake. Then around three months before the event, the situation changed. The fault region was going through a hot thermal anomaly state, breaking away from the previous cold anomaly state; i.e., the temperature of the fault zone underwent a change process from cold to hot status. After the earthquake, the value of the K tec in the fault region declined, exhibiting an increasingly discrete state in terms of spatial distribution. A similar phenomenon can be also observed in the Yutian earthquake event which occurred on the Altyn fault on 12 February 2014.
The above mentioned observed phenomenon appears to not be coincidental. In fact, we have investigated the characterization of TIR anomalies in tectonic activities before an earthquake event year via the TTIA method using four earthquake cases of magnitudes greater than Ms 7.0 over the Tibet Plateau since 2003, similar phenomena were observed in each; namely, before the large earthquake, the temperature of the surrounding rock underwent a heating to cooling process, while the temperature of the fault region underwent an inverse process. However, when we investigated earthquakes of magnitude Ms 6.0 or less, we did not observe this kind of very obvious regularity.
The possible explanation for the TIR anomalies in tectonic activities occurring before a large earthquake event is that it is the integration result of the green house air degassing from the Earth's crust and the stress of rocks. On one hand, green house air is able to warm up the air, thereby leading to three phase changes of water vapor, which can produce the variation of the air temperature. On the other hand, the compression and tension of rocks can result in LST variation through a series of complicated geophysical and geochemical processes.      Figure 17 shows the comparison between the MOD11A2 LST and air temperature observed by the ground site and the TIR anomalies value based on the TTIA method. It can be clearly seen that  Figure 17 shows the comparison between the MOD11A2 LST and air temperature observed by the ground site and the TIR anomalies value based on the TTIA method. It can be clearly seen that the MOD11A2 LST data are in high correspondence with the air temperature from the ground site, which can verify the efficiency of the MOD11A2 LST data. Meanwhile, the air temperature observed by the ground site is always higher than the MODIS LST, which can be attributed to the fact that the air temperature data are obtained during the day, while the MODIS LST used in the study are observed at night. The annual variation and very distinct seasonal characteristics of the original MOD11A2 LST data can be clearly observed, and this is the most important disturbance factor influencing the extraction of the tectonic TIR anomaly information. This is the reason for which we must first remove the influences of the solar radiation in the algorithm. Upon closer observation, we can see that the evolution process of the TIR anomalies based on the proposed TTIA algorithm is different from the original MODIS LST and ground air temperature, which indicates that the TTIA method can effectively remove the influence of solar radiation.

Comparison of the MOD11A2 LST and Ground Air Temperature Data, TIR Anomalies Value
Remote Sens. 2017, 9, x FOR PEER REVIEW 23 of 34 the MOD11A2 LST data are in high correspondence with the air temperature from the ground site, which can verify the efficiency of the MOD11A2 LST data. Meanwhile, the air temperature observed by the ground site is always higher than the MODIS LST, which can be attributed to the fact that the air temperature data are obtained during the day, while the MODIS LST used in the study are observed at night. The annual variation and very distinct seasonal characteristics of the original MOD11A2 LST data can be clearly observed, and this is the most important disturbance factor influencing the extraction of the tectonic TIR anomaly information. This is the reason for which we must first remove the influences of the solar radiation in the algorithm. Upon closer observation, we can see that the evolution process of the TIR anomalies based on the proposed TTIA algorithm is different from the original MODIS LST and ground air temperature, which indicates that the TTIA method can effectively remove the influence of solar radiation. Figure 17. Comparison of the temperature between from MODIS LST data and air temperature acquired by the Yutian meteorological ground station, and TIR anomaly based on the TTIA algorithm. The blue line denotes the daytime air temperature data observed by the ground site, the red line denotes the MODIS LST data at night, and the black line denotes the TIR anomaly based on the TTIA algorithm.

Comparison of the TTIA and RST Algorithms
The RST method is also an effective TIR anomalies extraction algorithm. It is based on a multitemporal analysis of the historical data set, which is defined as follows: where   ,  r x y represents the location coordinates on a satellite image;  t is the time of acquisition of the satellite image at hand, with t ∈ τ, where τ defines the homogeneous domain of the satellite imagery collected in the same time-slot (hour) of the day, and the same period (month) of the year;   ,   T r t denotes the difference between the current (   t t ) TIR signal value   ,  T r t at location r , and its spatial average    T t , computed in place on the image at hand, with the cloudy pixels discarded, and only the land pixels considered; ∆ ( ) and ∆ ( ) are the average and standard deviation values of ∆ ( , ), at location r , respectively, computed on cloud-free satellite records belonging to a selected homogeneous data-set (   t ). Although Equation (13) shows a certain similarity to Equation (14) in form, there were evident differences observed between the two. The input of the RST algorithm was the difference values between the LST and the mean value of all of the LST data in a spatial domain observed at the same time, the objective of which was to filter out the weather disturbances. However, the input of the TTIA algorithm was the , which was more closely related to the tectonic activities, due to the fact that the influences of solar radiation, atmosphere, and human activities had been filtered out. Figure 17. Comparison of the temperature between from MODIS LST data and air temperature acquired by the Yutian meteorological ground station, and TIR anomaly based on the TTIA algorithm. The blue line denotes the daytime air temperature data observed by the ground site, the red line denotes the MODIS LST data at night, and the black line denotes the TIR anomaly based on the TTIA algorithm.

Comparison of the TTIA and RST Algorithms
The RST method is also an effective TIR anomalies extraction algorithm. It is based on a multi-temporal analysis of the historical data set, which is defined as follows: where r ≡ (x, y) represents the location coordinates on a satellite image; t is the time of acquisition of the satellite image at hand, with t ∈ τ, where τ defines the homogeneous domain of the satellite imagery collected in the same time-slot (hour) of the day, and the same period (month) of the year; ∆T(r, t ) denotes the difference between the current (t = t ) TIR signal value T(r, t ) at location r, and its spatial average T(t ), computed in place on the image at hand, with the cloudy pixels discarded, and only the land pixels considered; µ ∆T(r) and σ ∆T(r) are the average and standard deviation values of ∆T(r, t), at location r, respectively, computed on cloud-free satellite records belonging to a selected homogeneous data-set (t ∈ τ).
Although Equation (13) shows a certain similarity to Equation (14) in form, there were evident differences observed between the two. The input of the RST algorithm was the difference values between the LST and the mean value of all of the LST data in a spatial domain observed at the same time, the objective of which was to filter out the weather disturbances. However, the input of the TTIA algorithm was the T tec , which was more closely related to the tectonic activities, due to the fact that the influences of solar radiation, atmosphere, and human activities had been filtered out.
The experimental results shown in Figure 18 exhibited the spatial distribution at Tibetan Plateau scope of the LST-based TIR anomalies with a contiguous group distribution characterization, rather than extending along the tectonic faults like the TTIA-based TIR anomalies shown in Figure 6. This was the main observed distinct difference between the TTIA and RST methods in regard to reflecting the spatial distribution of TIR anomalies on a large plateau scope. The extracted TIR anomalies space form prompted us to believe that the TTIA method is more sensitive to the tectonic TIR anomalies.
Remote Sens. 2017, 9, x FOR PEER REVIEW 24 of 34 The experimental results shown in Figure 18 exhibited the spatial distribution at Tibetan Plateau scope of the LST-based TIR anomalies with a contiguous group distribution characterization, rather than extending along the tectonic faults like the TTIA-based TIR anomalies shown in Figure 6. This was the main observed distinct difference between the TTIA and RST methods in regard to reflecting the spatial distribution of TIR anomalies on a large plateau scope. The extracted TIR anomalies space form prompted us to believe that the TTIA method is more sensitive to the tectonic TIR anomalies. Figure 18. TIR anomalies spatio-temporal presentation of the Ms 8.0 Wenchuan earthquake event based on the RST method on the Tibetan Plateau. Figure 19 shows the comparison of the spatio-temporal evolution of TIR anomalies extracted by the RST method and TTIA method respectively in the Longmenshan fault region six months to two months before the devastating Ms 8.0 Wenchuan earthquake event of 2008. From Figure 19, the following facts can be clearly observed. The spatial distribution of TIR anomalies based on the RST are roughly similar to those based on the TTIA method, such as in Figure 19a(3), b(3), yet nevertheless there was still some difference in the detailed section in the TIR anomalies between both methods, i.e., the distribution features of the clumps of TIR anomalies were more obvious in the RST method than the TTIA method. In addition, comparing Figure 19c(5-8) and d (5)(6)(7)(8), we can observe that, although the general distribution of the TIR anomalies was similar, the contrast between the positive TIR anomalies in the surrounding rocks region and negative TIR anomalies in the Longmenshan fault region based on the TTIA method was distinctly stronger than that based on the RST method, which demonstrated that the TTIA method was possibly more sensitive to the tectonic TIR anomalies signal.  Figure 19 shows the comparison of the spatio-temporal evolution of TIR anomalies extracted by the RST method and TTIA method respectively in the Longmenshan fault region six months to two months before the devastating Ms 8.0 Wenchuan earthquake event of 2008. From Figure 19, the following facts can be clearly observed. The spatial distribution of TIR anomalies based on the RST are roughly similar to those based on the TTIA method, such as in Figure 19a(3), b(3), yet nevertheless there was still some difference in the detailed section in the TIR anomalies between both methods, i.e., the distribution features of the clumps of TIR anomalies were more obvious in the RST method than the TTIA method. In addition, comparing Figure 19c(5-8),d (5)(6)(7)(8), we can observe that, although the general distribution of the TIR anomalies was similar, the contrast between the positive TIR anomalies in the surrounding rocks region and negative TIR anomalies in the Longmenshan fault region based on the TTIA method was distinctly stronger than that based on the RST method, which demonstrated that the TTIA method was possibly more sensitive to the tectonic TIR anomalies signal. In addition, it must be emphasized that the red dots indicate that earthquakes with magnitude between Ms6.0 to Ms6.5 took place in the corresponding period. The red arrow shows that an earthquake of magnitude greater than Ms6.5 occurred. The red solid arrow lines indicate that the earthquakes occurred inside the fault zone. The red dotted arrow lines indicate that the earthquake occurred outside the fault zone, yet the epicenter of the earthquake is close to the fault zone, so that they may have produced a certain impact on the tectonic thermal field in the fault zone.
Comparing Figures 20 and 21, we can observe some significant differences between the two methods: (1) The cluster based on the TTIA algorithm presents a more aggregate distribution trend than that based on the RST algorithm during an unperturbed period. This means that the TTIA algorithm has a lower false warning rate. For example, during the period of 2003 to 2004, when no large earthquake occurred, corresponding with the period marked by Figure 20a, in the region of Longmenshan Fault, all of the -mean is less than μ + 2σ based on the new algorithm. In contrast, during the same period, there are four occurrences where the value exceeded μ + 2σ, with the RST method, as shown in Figure 21, which indicates that the TTIA algorithm has a lower false alarm rate compared with the RST algorithm; (2) During the seismic periods, the warning rate of earthquakes of the TTIA algorithm is higher than that of the RST. For example, in the Longmenshan Fault zone, the Ms8.0 Wenchuan earthquake broke out on 12 May 2008 and the Ms7.0 Ya'an earthquake occurred on 20 April 2013. Before these two earthquakes, the value with the TTIA appeared to be more prominent than the RST, as exhibited in stages (b) and (c) of Figures 20 and 21. Specifically, during the period of one year before the Wenchuan earthquake, there were four instances where the value exceeded the value of μ + 2σ, based on the TTIA method. However, there were only three instances where value exceeded the value of μ + 2σ, as shown based on the RST method. Similarly, one year before the Ya'an earthquake, there were four instances where the value exceeded μ + 2σ, as shown by the TTIA method. In contrast, the value did not exceed μ + 2σ, as shown by the RST method. it must be emphasized that the red dots indicate that earthquakes with magnitude between Ms 6.0 to Ms 6.5 took place in the corresponding period. The red arrow shows that an earthquake of magnitude greater than Ms 6.5 occurred. The red solid arrow lines indicate that the earthquakes occurred inside the fault zone. The red dotted arrow lines indicate that the earthquake occurred outside the fault zone, yet the epicenter of the earthquake is close to the fault zone, so that they may have produced a certain impact on the tectonic thermal field in the fault zone. Comparing Figures 20 and 21, we can observe some significant differences between the two methods: (1) The cluster based on the TTIA algorithm presents a more aggregate distribution trend than that based on the RST algorithm during an unperturbed period. This means that the TTIA algorithm has a lower false warning rate. For example, during the period of 2003 to 2004, when no large earthquake occurred, corresponding with the period marked by Figure 20a, in the region of Longmenshan Fault, all of the K tec -mean is less than µ + 2σ based on the new algorithm. In contrast, during the same period, there are four occurrences where the K tec value exceeded µ + 2σ, with the RST method, as shown in Figure 21, which indicates that the TTIA algorithm has a lower false alarm rate compared with the RST algorithm; (2) During the seismic periods, the warning rate of earthquakes of the TTIA algorithm is higher than that of the RST. For example, in the Longmenshan Fault zone, the Ms 8.0 Wenchuan earthquake broke out on 12 May 2008 and the Ms 7.0 Ya'an earthquake occurred on 20 April 2013. Before these two earthquakes, the K tec value with the TTIA appeared to be more prominent than the RST, as exhibited in stages (b) and (c) of Figures 20 and 21. Specifically, during the period of one year before the Wenchuan earthquake, there were four instances where the K tec value exceeded the value of µ + 2σ, based on the TTIA method. However, there were only three instances where K tec value exceeded the value of µ + 2σ, as shown based on the RST method. Similarly, one year before the Ya'an earthquake, there were four instances where the K tec value exceeded µ + 2σ, as shown by the TTIA method. In contrast, the K tec value did not exceed µ + 2σ, as shown by the RST method.   Figure 22(a1),(b1), it can be seen that original landsurface temperature is mainly controlled by the topography because there is a good correspondence between the terrain of valley and high landsurface temperature. Figure 22(a2), (b2) also shows a similar pattern. This indicates that the topography is the important leading factor of the original landsurface temperature. Figure 22(c1), (c2) presented the TIR anomalies extracted by the TTIA algorithm. It can be obviously seen that there is no correspondence between the TTIA-based TIR anomalies and the terrain, which demonstrates that the TTIA algorithm is able to remove the influence of the terrain, and simultaneously highlight the active tectonic lineations characteristic. Particularly, comparing Figure 22(a2),(c2), it can be easily seen that there are very distinct difference between the value of the DEM and TIR anomalies in the Longmenshan fault zone, which indicates once more that the TTIA algorithm is able to remove the influence of the topographic effect of ridges/valley in such a mountainous region and reflect the tectonic features well.   Figure 22(a1),(b1), it can be seen that original landsurface temperature is mainly controlled by the topography because there is a good correspondence between the terrain of valley and high landsurface temperature. Figure 22(a2), (b2) also shows a similar pattern. This indicates that the topography is the important leading factor of the original landsurface temperature. Figure 22(c1), (c2) presented the TIR anomalies extracted by the TTIA algorithm. It can be obviously seen that there is no correspondence between the TTIA-based TIR anomalies and the terrain, which demonstrates that the TTIA algorithm is able to remove the influence of the terrain, and simultaneously highlight the active tectonic lineations characteristic. Particularly, comparing Figure 22(a2),(c2), it can be easily seen that there are very distinct difference between the value of the DEM and TIR anomalies in the Longmenshan fault zone, which indicates once more that the TTIA algorithm is able to remove the influence of the topographic effect of ridges/valley in such a mountainous region and reflect the tectonic features well.  Figure 22 shows the DEM (Digital Elevation Model) data, original MODIS LST data and the TIR anomalies based on TTIA algorithm. Comparing the Figure 22(a1),(b1), it can be seen that original landsurface temperature is mainly controlled by the topography because there is a good correspondence between the terrain of valley and high landsurface temperature. Figure 22(a2),(b2) also shows a similar pattern. This indicates that the topography is the important leading factor of the original landsurface temperature. Figure 22(c1),(c2) presented the TIR anomalies extracted by the TTIA algorithm. It can be obviously seen that there is no correspondence between the TTIA-based TIR anomalies and the terrain, which demonstrates that the TTIA algorithm is able to remove the influence of the terrain, and simultaneously highlight the active tectonic lineations characteristic. Particularly, comparing Figure 22(a2),(c2), it can be easily seen that there are very distinct difference between the value of the DEM and TIR anomalies in the Longmenshan fault zone, which indicates once more that the TTIA algorithm is able to remove the influence of the topographic effect of ridges/valley in such a mountainous region and reflect the tectonic features well.   Figure 1 with F2) is a active tectonic zone. Meanwhile the TIR anomalies based on the TTIA algorithm present a very obvious distribution characteristic along the active tectonic zone no matter the event year or the undisturbance year. Previous studies demonstrate that the TIR anomalies correspond with main tectonic lines [2,8,21,43], due to the fact that the active tectonic lines made up of the faults are often the channel of greenhouse gases, such as CO2, NH4 and so on, degassing from rocks under stress, and as well as the participation of ground water as a possible cause for generation of TIR anomalies or p-hole activation in the stressed rock volume and their further recombination at the rock-air interface [9,23,[50][51][52][53][54][55][56]. The study results demonstrate that the TTIA algorithm is an effective algorithm to extract the TIR anomalies ralated to the tectonic zone. However, further comparing the Figure 23(C1),(D1), we can find that in an event year, the TIR anomalies are more dramatic than that in an undisturbance year.   Figure 1 with F2) is a active tectonic zone. Meanwhile the TIR anomalies based on the TTIA algorithm present a very obvious distribution characteristic along the active tectonic zone no matter the event year or the undisturbance year. Previous studies demonstrate that the TIR anomalies correspond with main tectonic lines [2,8,21,43], due to the fact that the active tectonic lines made up of the faults are often the channel of greenhouse gases, such as CO 2 , NH 4 and so on, degassing from rocks under stress, and as well as the participation of ground water as a possible cause for generation of TIR anomalies or p-hole activation in the stressed rock volume and their further recombination at the rock-air interface [9,23,[50][51][52][53][54][55][56]. The study results demonstrate that the TTIA algorithm is an effective algorithm to extract the TIR anomalies ralated to the tectonic zone. However, further comparing the Figure 23(C1),(D1), we can find that in an event year, the TIR anomalies are more dramatic than that in an undisturbance year. In a word, through the comparison of the TTIA algorithm-based TIR anomalies and original MODIS LST data, DEM data, GPS displacement data and as well as the active tectonic lines extension direction, we can know that the TTIA algorithm can reduce the disturbance of the atectonic factors, such as solar radiation, air circulation, human activities, and terrain. The study result shows that the TTIA algorithm is effective for the characterization of thermal infrared anomalies in tectonic activities. Moreover, the dangerous region which has high possibility to meet an impending earthquake can be sketched by calculating the gentroid of the TIR anomalies based on this TTIA algorithm. It should be noted that the outline of earthquake danger zone based only on TIR anomalies is not enough. Other pre-earthquake anomalies such as the electromagnetic anomalies in the atmospheric ionosphere, extracted by the Swarm satellite [93] and DEMETER satellite [94] should been incorporated.

Conclusions
In this study, a novel TTIA (tectonic thermal infrared anomalies) method was proposed to characterize earthquake TIR anomalies. The Tibetan Plateau was chosen as the study area. Then, using the MODIS land surface temperature (LST) products MOD11A2 via the proposed method, this study observed the spatio-temporal evolution of the TIR anomalies over the plateau from 2003 to 2015. A method that combined the typical earthquake cases, and analyzed the TIR anomalies signals of the pre-seismic, co-seismic, and post-seismic periods, was implemented, which allowed for a deeper understanding of the earthquakes' TIR anomalies features. The TIR anomalies rule at a plateau scale and fault zone scale was examined prior to the onset of the major earthquake events in In a word, through the comparison of the TTIA algorithm-based TIR anomalies and original MODIS LST data, DEM data, GPS displacement data and as well as the active tectonic lines extension direction, we can know that the TTIA algorithm can reduce the disturbance of the atectonic factors, such as solar radiation, air circulation, human activities, and terrain. The study result shows that the TTIA algorithm is effective for the characterization of thermal infrared anomalies in tectonic activities. Moreover, the dangerous region which has high possibility to meet an impending earthquake can be sketched by calculating the gentroid of the TIR anomalies based on this TTIA algorithm. It should be noted that the outline of earthquake danger zone based only on TIR anomalies is not enough. Other pre-earthquake anomalies such as the electromagnetic anomalies in the atmospheric ionosphere, extracted by the Swarm satellite [93] and DEMETER satellite [94] should been incorporated.

Conclusions
In this study, a novel TTIA (tectonic thermal infrared anomalies) method was proposed to characterize earthquake TIR anomalies. The Tibetan Plateau was chosen as the study area. Then, using the MODIS land surface temperature (LST) products MOD11A2 via the proposed method, this study observed the spatio-temporal evolution of the TIR anomalies over the plateau from 2003 to 2015. A method that combined the typical earthquake cases, and analyzed the TIR anomalies signals of the pre-seismic, co-seismic, and post-seismic periods, was implemented, which allowed for a deeper understanding of the earthquakes' TIR anomalies features. The TIR anomalies rule at a plateau scale and fault zone scale was examined prior to the onset of the major earthquake events in the area. Then, after analyzing this study's experimental results, the following conclusions were drawn: 1.
The obtained TIR anomalies based on the new algorithm showed an obvious spatial distribution characteristic along the main faults on the plateau. Therefore, it can be proved that the proposed algorithm had distinctive advantages in removing or weakening the disturbances of the atectonic factors, and was therefore very effective in extracting the tectonic TIR anomalies signals.

2.
The seismogenic zone was found to be a more effective observation scope for the deeper understanding of the mid-and short-term seismogenic and crust stress change processes. 3.
The movement trace of the centroid of the TIR anomalies over the entire plateau was helpful in judging the approximate dangerous tectonic regions where major earthquakes may occur in the future.

4.
At the observe scale of earthquake generating fault zone, before the great earthquake, the fluctuations of the K tec value are significantly more volatile than those in an aseismic period, which indicates that the TIR anomalies in tectonic activities in the event year are more active than those in non-event years.