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

Abstract

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 deficiencies of the classic algorithm, which is currently used for noise filtering 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.

1. 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,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43]. 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 ²²² 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 weather interferences on the extractions of the TIR anomalies [60]; (2) The algorithms, which are based on background analyses [61,62,63,64], such as robust satellite techniques (RST) [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64], are generally based on the multi-temporal analyses of the historical data sets acquired by satellite observations under similar observational conditions (for example, the same month of the year; same hour of the day; same sensor, and so on). In these methods, the spatial-temporal anomalies are always identified using a comparison of the preliminarily computed signal behaviors, such as background information; (3) The algorithms are based on signal analyses [65,66,67,68,69,70,71,72], such as a wavelet decomposition method [65]; wavelet packet method [66]; and power spectrum method [69,70,71,72]. These methods extract the useful characteristic signals in a frequency domain.

Currently, many arguments still exist in regard to the efficiency of algorithms’ abilities to indicate tectonic activities. Although these methods are known to be able to remove the impacts of solar radiation, they cannot effectively address the filtering of atmospheric disturbances and human activities. Therefore, in order to extract as much possible TIR anomalies information, which is highly related to tectonic activities, this study proposed a novel algorithm to characterize tectonic TIR anomalies, which is referred to in this study as TTIA (tectonic thermal infrared anomalies). Two important devastating earthquakes (Ms 8.0 Wenchuan earthquake in 2008, and Ms 7.3 Yutian earthquake in 2014) were selected as the study cases. The spatial-temporal evolution of the tectonic TIR anomalies based on the proposed TTIA algorithm during the entire time period (2003 to 2015) were presented. In this study, the tectonic TIR anomalies from the macroscopic scope were successfully observed, and some new findings were obtained regarding the TIR anomalies evolution characteristics prior to major earthquake events. The results of this research study may potentially provide important references for predicting the dangerous regions where significant earthquakes may occur in the future.

The main objectives of this study were as follows: (1) To propose a novel algorithm to characterize the TIR anomalies which are highly related to tectonic activities; and (2) To demonstrate and summarize the tectonic TIR anomalies evolution rule of two previous major earthquakes, and provide ideas for outlining potentially dangerous regions where major earthquakes may occur in the future.

The remainder of this research document is organized as follows: In Section 2, the MODIS LST data used in this study, along with the methodology, are described. In Section 3, the experimental results are presented. The discussion and conclusions reached in this study are detailed in Section 4 and Section 5, respectively.

2. Methodology

2.1. Study Area and the Two Examined Earthquake Cases

The Qinghai-Tibet Plateau was taken as the study area. As a whole, the Qinghai-Tibet Plateau was actually formed by isostatic mountain building processes beginning in 3.6 Ma, which were accompanied by the isostatic subsidence of the peripheral basins and depressions, and major changes in the topography and environment of the area. A series of pulsating mountain building events occurred at 3.6 Ma, 2.5 Ma, 1.8 to 1.2 Ma, 0.9–0.8 Ma, and 0.15–0.12 Ma, which were associated with the equilibrium of gravity and the isostactic adjustments of the crustal materials. There events led to the formation of a composite orogenic belt in which the originally relatively independent Himalayas, Gangdise, Tanghla, Longmenshan, Kunlun, Altyn Tagh, and Qilian Mountains became unified. The result was the formation of the complete Qinghai-Tibet Plateau with a unified mountain root following the Miocene uplift of the entire plateau [73]. With its very thick continental crust (70 km thickness on average; double the normal thickness of the crust), and the dramatic uplift of the late Pleistoeene period (now with an average elevation of between 4500 and 5000 m), along with its sustained and strong tectonic deformations [74], the Qinghai-Tibet Plateau has become one of the most unique geological areas in the world. It spans across Eurasia, and is an important part of the Alps-Himalayan tectonic belt. The stratum development of the plateau is now complete. The drastic lifting movement is known as the “the roof of the world”, where a series of large-scale and active fault zones have been formed since the Cenozoic. These fault zones include the following: the Himalaya main thrust belt (F1); Karakorum-Jiali fault zone (F2); Mani-Yushu-Xianshuihe fault zone (F3); Kunlun-Maqin fault zone (F4); Altyn-Haiyuan fault zone (F5); and the Jinshajiang-Honghe fault zone (F6), from the south to the north. These six main fault zones divide the Tibetan Plateau into six different active blocks as follows: The Chuandian (B1); Lasa (B2); Qiangtang (B3); Bayankala (B4); Chaidamu (B5); and Qilian (B6) [69], as shown in Figure 1.

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 Ms 7.0 consisted of the Ms 8.0 Wenchuan earthquake of 2008; the Ms 7.3 Yushu earthquake of 2010; the Ms 7.0 Yaan earthquake of 2013; and the Ms 7.3 Yutian earthquake which occurred in 2014. These earthquake events caused enormous economic losses and casualties.

On 12 May 2008, at 14:28 local time (06:28 UTC), the devastating Ms 8.0 Wenchuan (Sichuan, China) earthquake ruptured the Longmen Shan fault zone on the eastern edge of China’s Tibetan Plateau (illustrated in the yellow square box labelled 2 in 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.

2.2. 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. Data set of the MOD11A2.

Name	Data Type	Effective Numerical Range	Unit	Filling Values	Calibration Coefficient
LST_Day_1 km: 8-Day daytime 1 km grid LST	16-bit unsigned int	7500–65535	K	0	0.02
QC_Day: Quality control for daytime LST and emissivity	8-bit unsigned int	0–255
Day_view_time: Average time of daytime LST observation	8-bit unsigned int	0–240	h	255	0.1
Day_view_angle: Average view zenith angle of daytime LST	8-bit unsigned int	0–130	Degree	255	1(-65)
LST_Day_1 km:8-Day nighttime 1 km grid LST	16-bit unsigned int	7500–65535	K	0	0.02
QC_Day: Quality control for nighttime LST and emissivity	8-bit unsigned int	0–255
Day_view_time: Average time of nighttime LST observation	8-bit unsigned int	0–240	h	255	0.1
Day_view_angle: Average view zenith angle of nighttime LST	8-bit unsigned int	0–130	Degree	255	1 (−65)
Emis_32: Band32 emissivity	8-bit unsigned int	1–255		0	0.002 (+0.49)
Emis_31: Band31 emissivity	8-bit unsigned int	1–255		0	0.002 (+0.49)
Clear_sky_days: the days in clear sky conditions and with valid LST	8-bit unsigned int	1–255		0
Clear_sky_nights: the nights in clear sky conditions and with valid LST	8-bit unsigned int	1–255		0

. 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.

2.3. 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 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 ③ 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 ⑤ 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.

2.3.1. 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 $\rho \left(x\right)\ge 0$. Where $\rho \left(x\right)$ 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\left[a,b\right]$, then the norm of $f(x)$ can be obtained using Equation (1) [82] as follows:

$${\Vert f(x)\Vert }^2={\displaystyle {\int }_a^b\rho }\left(x\right)f^2\left(x\right)dx$$

(1)

The nature of the optimal approximation in $L_2$ is that, for any $f\left(x\right)\in \left(a,b\right)$, function $\phi \left(x\right)$ should be looked for on the basis space $M=\left\{{\phi }_1,{\phi }_2,\dots ,{\phi }_n\right\}$ as follows:

$$\phi \left(x\right)={\displaystyle {\int }_{i=1}^n{c}_i{\phi }_i}$$

(2)

which is subject to the following:

$${\displaystyle {\int }_a^b\rho \left(x\right){[f\left(x\right)-\phi \left(x\right)]}^{2}dx\to \min }$$

(3)

Equation (3) indicated that the square of the error function reached a minimum value after solving the integral. Therefore, $\phi \left(x\right)$ 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 $\phi \left(x\right)$ will be the optimal Fourier approximation of $\phi \left(x\right)$.

➢ 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:

$$y={\displaystyle \sum _{i=1}^m\left[a_i\sin \left({\omega }_{i}t+b_i\right)+c_i\cos \left({\omega }_{i}t+b_i\right)\right]}+\epsilon$$

(4)

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; ${\omega }_i$ is the phase of the $i_{th}$ harmonic; and $\epsilon$ 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:

$${\displaystyle \sum {\epsilon }^2={\displaystyle \sum _{j=1}^n{\left\{y_j-{\displaystyle \sum _{i=1}^m\left[a_i\sin \left({\omega }_{i}t+b_i\right)+c_i\cos \left({\omega }_{i}t+b_i\right)\right]}\right\}}^2}}$$

(5)

2.3.2. 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.

➢ Theory of the 2-D wavelet transform

● 2-D continuous wavelet transform

The 2-D continuous wavelet transform was defined as follows [7]:

$$W_f\left(a,b_x,b_y\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}f\left(x,y\right){\mathsf{\Psi }}_{a,b_x,b_y}\left(x,y\right)dxdy}}$$

(6)

$${\mathsf{\Psi }}_{a,b_x,b_y}\left(x,y\right)=\frac{1}{a}\mathsf{\Psi }\left(\frac{x-b_x}{a},\frac{y-b_y}{a}\right)$$

(7)

Its inverse was given by Equation (8):

$$f\left(x,y\right)=\frac{1}{C_{\mathsf{\Psi }}}{\displaystyle \underset{0}{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{W}_f\left(a,b_x,b_y\right){\mathsf{\Psi }}_{a,b_x,b_y}\left(x,y\right)d{b}_{x}d{b}_y\frac{da}{a^3}}}}$$

(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$, $\mathsf{\Psi }\left(x,y\right)$ was elongated; and when $a>1$, $\mathsf{\Psi }\left(x,y\right)$ 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]:

$${\mathsf{\Psi }}_{m,n}\left(t\right)=\frac{1}{\sqrt{a_0^m}}\mathsf{\Psi }\left(\frac{t-n{b}_0{a}_0^m}{a_0^m}\right)=a_0^{-\frac{m}{2}}\mathsf{\Psi }\left(a_0^{-m}t-n{b}_0\right)$$

(9)

where $m,n$ are the integer,$a_0\ne 1$ and $b_0\ne 0$. The corresponding discrete wavelet transform can then be defined as follows:

$$\langle f,{\mathsf{\Psi }}_{m,n}\rangle =a_0^{-\frac{m}{2}}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}f\left(t\right)\mathsf{\Psi }\left(a_0^{-m}t-n{b}_0\right)dt}$$

(10)

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):

$$T_{tec}=T_{L4}-T_{L9}$$

(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.

2.3.3. Expressions of the Tectonic Thermal Infrared Anomalies

➢ “$\mathit{k}\mathit{\sigma }$”rule and offset index ${\mathit{K}}_{\mathit{t}\mathit{e}\mathit{c}}$

In this study, the “$k\sigma$”rule was used to determine the signals of the thermal anomalies. The “$k\sigma$” principle was proposed initially by Wright in 1884 [90,91,92,93]. This simple and basic principle is as follows:

$$\left|x_i-\overline{x}\right|\ge k\sigma$$

(12)

where x and σ are the mean and standard deviation of the result x, respectively. Here the “$k\sigma$” 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 $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 ${\sigma }_{T_{tec}}\left(r_i\right)$.

$$K_{tec}=\frac{T_{tec}\left(r_i,t\right)-{\mu }_{T_{tec}}\left(r_i\right)}{{\sigma }_{T_{tec}}\left(r_i\right)}$$

(13)

where $r_i$ denotes the location of pixel $i$ on the remote sensing image, with the abscissa and ordinate $\left(x_i,y_i\right)$; $t$ represents the time of the acquisition of the satellite image, with $t\in \tau$, where $\tau$ 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 ${\mu }_{T_{tec}}$ and ${\sigma }_{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 ${\mu }_{T_{tec}}$, $\mathrm{and}{\sigma }_{T_{tec}}$ in the Equation (13) has distinct $T_{tec}$ basis, compared with the $\mu$ $\mathrm{and}\sigma$ in the RST method.

3. 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}\ge \mu +\sigma$. In this study, $K_{mean}$ is the spatial mean of the $K_{tec}$ value image in a period; $\mu$ is the mean of the $K_{tec}$ value images in the multi-periods; and $\sigma$ 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(6)–(8).

Figure 8 shows the enlarged images from Figure 7. The TIR anomalies were derived by means of the proposed algorithm in the Lonmenshan fault zone where the Ms 8.0 Wenchuan earthquake occurred. It can be seen that the TIR anomalies based on the TTIA method are in high correspondence with the tectonic line.

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.

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.

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 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 spatio-temporal 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 Figure 13, Figure 14, Figure 15 and Figure 16. 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.

4. Discussion

4.1. Comparison of the MOD11A2 LST and Ground Air Temperature Data, TIR Anomalies Value

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.

4.2. 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:

$${\otimes }_{\Delta T}\left(r,t^{\prime }\right)=\frac{\Delta T\left(r,t^{\prime }\right)-{\mu }_{\Delta T}\left(r\right)}{{\sigma }_{\Delta T}\left(r\right)}$$

(14)

where $r\equiv \left(x,y\right)$ represents the location coordinates on a satellite image; $t^{\prime }$ 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; $\Delta T\left(r,t^{\prime }\right)$ denotes the difference between the current ($t=t^{\prime }$) TIR signal value $T\left(r,t^{\prime }\right)$ at location $r$, and its spatial average $T\left(t^{\prime }\right)$, computed in place on the image at hand, with the cloudy pixels discarded, and only the land pixels considered; ${\mu }_{\Delta T\left(r\right)}$ and ${\sigma }_{\Delta T\left(r\right)}$ are the average and standard deviation values of $\Delta T\left(r,t\right)$, at location $r$, respectively, computed on cloud-free satellite records belonging to a selected homogeneous data-set ($t^{\prime }\in \tau$).

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.

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–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.

The calculated results of TTIA and RST in the Longmenshan Fault zone are presented in a time series, shown as Figure 20 and Figure 21. Each black or red dot indicates the mean value of TIR anomalies throughout the fault zone (${10}^{°}\times {10}^{°}$) in a certain 8-day period between 2003 and 2015. In addition, 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 Figure 20 and Figure 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 Figure 20 and Figure 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.

4.3. Comparison of TIR Anomalies Based on TTIA and Tectonic Lineations, Topographic Effect

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 23 shows the DEM data, original MODIS LST data, TIR anomalies based on the TTIA algorithm in the event year and the undisturbance year (no Ms ≥ 5 earthquakes), GPS displacement data, and tectonic fault vector data over the Tibet-plateau. Comparing Figure 23(A1),(B1), and Figure 23(A2),(B2), same rule can be obtained; i.e., the original landsurface temperature on the whole Tibet plateau is controlled by the topography. The black arrows in Figure 23(C1) indicate the vector of the GPS displacement. Observing the GPS displacement, we can know that the Karakorum-Jiali fault zone (shown in 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₂, NH₄ 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.

5. 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:

• 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.
• 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.
• 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.
• 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.