Global land surface temperature influenced by vegetation cover and PM2.5 from 2001 to 2016 Global land surface temperature influenced by vegetation cover and PM2.5 from 2001 to 2016. Remote

: Land surface temperature (LST) is an important parameter to evaluate environmental changes. In this paper, time series analysis was conducted to estimate the interannual variations in global LST from 2001 to 2016 based on moderate resolution imaging spectroradiometer (MODIS) LST, and normalized difference vegetation index (NDVI) products and ﬁne particulate matter (PM 2.5 ) data from the Atmospheric Composition Analysis Group. The results showed that LST, seasonally integrated normalized difference vegetation index (SINDVI), and PM 2.5 increased by 0.17 K, 0.04, and 1.02 µ g/m 3 in the period of 2001–2016, respectively. During the past 16 years, LST showed an increasing trend in most areas, with two peaks of 1.58 K and 1.85 K at 72 ◦ N and 48 ◦ S, respectively. Arctic. relations in the Southern Hemisphere, negative correlations were detected. The impact of PM 2.5 on LST was more complex. On the whole, LST increased with a small increase in PM 2.5 concentrations but decreased with a marked increase in PM 2.5 . The study provides insights on the complex relationship between vegetation cover, air pollution, and land surface temperature.

to the growing emissions of greenhouse gas [45]. Previous studies indicated that air pollution might produce a fog-like haze which could reduce the solar radiation and might affect the precipitation and microphysical properties of clouds [42,46,47]. For example, haze caused a 0.7 K temperature difference between urban and rural area at nights in semi-arid cities of China [48]. However, most of current studies about PM 2.5 researched the influence on meteorological conditions at a regional scale. Only a few studies have researched the spatio-temporal characteristics of PM 2.5 at a global scale.
Ground-based measurements could provide the most reliable values of land surface temperatures. However, due to limitations caused by sparsity and the inhomogeneous distribution of observation stations, these measurements are insufficient for global-scale research. For example, in some regions, such as Africa and the Antarctic, where very few stations are available and therefore this will prevent us from understanding the changing surface temperatures and extreme weather [17]. New observation methods are still needed. With the development of remote sensing technology, remote sensing images have been increasingly used on assessing changes in large spatiotemporal scales. The modern satellite techniques can provide continuous and temporally repetitive images and long-term observations at the global scale, and it provides the only possibility to quickly measure really information of the entire globe [7,24,49]. The moderate resolution imaging spectroradiometer (MODIS) is an important sensor that is installed on the Terra and Aqua satellites of the Earth Observing System (EOS) and operated by the National Aeronautics and Space Administration (NASA). Due to combination of the global coverage, intermediate spatial resolution, higher temporal resolution, and abundant products, MODIS data have been widely used in environmental research [50][51][52].
In short, previous studies on LST mainly focused on retrieval algorithm, accuracy verification and urban heat island at the regional scale. There is a lack of comprehensive understanding among LST, vegetation cover and PM 2.5 at a global scale. Therefore, it is necessary to explore the influence vegetation cover and PM 2.5 on LST. In this study, Terra MODIS LST and NDVI data were obtained from the Land Processes Distributed Active Archive Center (LP DAAC) managed by the NASA Earth Science Data and Information System (ESDIS) project. PM 2.5 data were provided by Dalhousie University Atmospheric Composition Analysis Group. Multitemporal data were used to analyze the spatiotemporal variations of LST. The effects of vegetation cover and PM 2.5 on LST were estimated during the period of 2001-2016 at a global scale. The objectives of this study are (1) to determine the spatiotemporal changes in LST, vegetation cover and PM 2.5 concentrations; (2) to explore the relationship between LST and vegetation cover changes measured by the seasonally integrated normalized difference vegetation index (SINDVI) and atmospheric contamination variations measured by PM 2.5 concentrations; and (3) to analyze the factors influencing LST.

Data and Processing
2.1.1. LST LST data were collected by the MODIS sensor, a key instrument of the Terra and Aqua satellites. The LST products were extracted from the MOD11C3 collection 6 dataset from 1 January 2001 to 31 December 2016, which provided monthly composited and averaged LST at a 0.05 • spatial resolution (5.6 km × 5.6 km at the Equator), and obtained from Earthdata Search (https://search.earthdata.nasa. gov/). There are 17 scientific data sets (SDS) in this product. The band named LST_Day_CMG was extracted as day time LST to analyze the variation of LST. The standard deviation of the MODIS V006 LST errors were within 1 K [53]. The mean value composite (MVC) method was used to calculate annual average values of LST in the period of 2001-2016.

NDVI
The Terra MODIS NDVI products were extracted from MOD13C2 collection 6 datasets, which were composed of 16-day global composite images from January 2001 to December 2016, and provided by the Earthdata Search, with a spatial resolution at 0.05 • . It is derived from the red and near-infrared reflectance ratio as follows [23]: where NIR and R are the amounts of near-infrared (the reflectance of band 2 of MODIS) and red light (the reflectance of band 1 of MODIS).
Although the product has eliminated some clouds and bidirectional reflectance of ground objects [54], MODIS NDVI time series data are still influenced by the clouds and atmospheric water vapor. A simple but efficient method, the mean-value iteration filter, was used to reduce the noise and reconstruct a high quality of NDVI time series by using the following equation [55]: where NDVI i represents the ith monthly NDVI, with range from 1 to 192 (16 years). A small percentage of the multiyear NDVI average of each pixel was set at a threshold (∆). The results of several test showed that 10% was the optimal threshold in this study. When ∆ i > ∆, the NDVI i was replaced by (NDVI i−1 + NDVI i+1 )/2. The procedure is iterated until all ∆ i are within the given threshold. The seasonal characteristics of plants can cause substantial differences in vegetation cover between seasons. Thus, the seasonally integrated normalized difference vegetation index (SINDVI), sum of NDVI values exceeding the threshold (commonly defined as NDVI > 0.1) [56][57][58], was used to monitor the interannual change of vegetation. It can identify vegetation distribution and growth state simultaneously. Due to the large influence of intermittent climatic factors such as precipitation on vegetation growth, compared with NDVI, SINDVI could better describe the interannual variation [59].
2.1.3. PM 2.5 PM 2.5 data were significantly consistent (R 2 = 0.81) with the monitoring results of PM 2.5 concentrations [60] obtained from the Atmospheric Composition Analysis Group (http://fizz.phys. dal.ca/~{}atmos/martin/) with a spatial resolution of 0.1 • . To match the resolutions of the LST and SINDVI datasets, the data were resampled to 0.05 • latitude × 0.05 • longitude using the nearest neighbor method. The unit for these data was µg/m 3 . These data only included the region spanning from 69.8 • N to 54.8 • S and 180 • E to 180 • W. The objective of this study was to detect global change trends and impact of PM 2.5 on LST. Because of lower PM 2.5 concentrations in the Arctic and Antarctica, these data could meet the target of this study.

Inter-Annual Variation Analysis
To determine the spatiotemporal differences, the annual change tendency was quantified at a global scale. Monthly time series data were processed to form yearly time series data; in this way, for each pixel the time series has a length of 16, spanning from 2001 to 2016, and the number of pixels in each band are 3600 × 7200. Linear regression and ordinary least squares methods were used to estimate the rate and range of changes at regional and global scales [9,57,61,62]. Random noises within 1K were added to LST data to test the sensitivity of constant biases. It was worth nothing that the slope and range were not affected by possible constant biases present in the data. Therefore, the slope and range were used to assess the interannual variations in LST, SINDVI, and PM 2.5 . The slope represents the change rate of each pixel from 2001 to 2016 calculated by the ordinary least-squares estimation via linear regression as follows: where n represents the length of the time series (in this research, n = 16); A i represents the ith year for SINDVI, PM 2.5 or the average LST; and the value of i ranges from 1 to n. The range can evaluate the magnitude of change. Each pixel value represents the total variation during the period of 2001-2016.
where n is similar to the definition of slope and represents the length of the time series. If slope > 0, it means that there is an increasing trend during the past 16 years, whereas the opposite means that it is decreasing. A large absolute value indicates a dramatic change, while the absolute value close to 0 indicates little variation.

Correlation Analysis
For clear descriptions of the relationships between LST with SINDVI and PM 2.5 , Pearson's correlation coefficients (r) was used to indicate the correlation. It was conducted as follows: where x i and y i represent LST and either SINDVI or PM 2.5 in the ith year, respectively; and x and y represent the corresponding averages over 16 years. At the same time, partial correlations were conducted to assess the dominant factor of each pixel. The partial correlation analysis refers to the process of removing the influence of the third variable when two variables are related to the third variable at the same time, and it only analyzes the correlation between the other two variables [56]. When the degree of correlation between LST and SINDVI was estimated, the influence of PM 2.5 could be removed. The partial correlation coefficient was calculated as follows: r xy(z) = r xy − r xz r yz where r xy(z) is the partial correlation coefficient, representing the degree of correlation between x and y without z; and r xy , r xz and r yz are simple correlation coefficients. The C program was used to calculate the average LST for each year from 192 scene images, reconstruct the NDVI, summarize the NDVI during growing season in order to obtain the annual SINDVI and analyze the correlation. The Mollweide projection, an equal-area pseudo-cylindrical map projection, was used to count regional areas.

Spatiotemporal Changes in LST, SINDVI, and PM 2.5 Concentrations
The global land surface temperature presented an increasing trend, with a global mean rise of 0.17 K from 2001 to 2016 (Table 1). LST increased at various degrees across 62.4% of the land area, with the exception in some regions such as the Antarctic, central North America, Central Asia, and eastern Australia (Figure 1a). An LST increase of 15 K was found in the Aral Sea, where water level decreased by 23 m, surface area decreased by 74%, and lake volume decreased by 90% [63], experiencing extreme desiccation [62]. The temperature of most regions increased within the range between 0 and 3 K, which accounts for 61.9% of the global land area. The Arctic experienced a marked warming trend, especially in northern Russia; however, the Antarctic, an area of approximately 14 million km 2 , showed a cooling trend. In the mid-latitudes of the Northern Hemisphere, such as southeast Canada, northern United States, eastern Kazakhstan, Mongolia, and Northeast China, dramatic decreases in LST were detected.
Variable solar radiation causes different climate zones and LST along the latitude. As shown in Figure 1a and Table 1, LST had an increasing trend, except for the areas at high latitudes in the Southern Hemisphere. In the high latitudes of the Northern Hemisphere, the change range increased firstly before a sharp decline. LST increased by an average of 0.30 K in the mid-latitudes of northern hemisphere, and it fluctuated slightly in the low-latitude areas. Compared with the Northern Hemisphere, there were two distinct peaks both at low-and mid-latitudes, and the larger one was an increase of 1.85 K at 48 • S.
SINDVI presented a greening trend all over the world [57], with an average value of 0.04 during the period of 2000-2016. Approximately 58.4% of the global area showed an increasing trend. Arctic, Antarctic, and high mountains were coved by snow for most time of the year. Approximately 9.8% of the global land was covered with snow all the year round, and SINDVI had no change in these regions. About 31.8% of other areas experienced a downward trend. Approximately 54.5% of the area's increase range of SINDVI was between 0 and 1, and 40.5% ranged from −1 to 0. The vegetation cover in central and southeastern China, northwestern India, and eastern Australia increased markedly. In addition, an increase appeared in Europe and most parts of North America. However, the SINDVI decreased significantly in eastern Brazil, northern Argentina, northern Kazakhstan, southeastern Africa, and western Australia.
The latitudinal SINDVI varied from −0.11 to 0.33, with an obvious increasing trend in the Northern Hemisphere, which indicated the great increase in the vegetation coverage in the Northern Hemisphere ( Figure 1b and Table 1). SINDVI in Southern Hemisphere (between 15 • S to 30 • S) also increased, with fluctuation within ±0.15. However, a dramatic decrease (from 0.32 to −0.09) was found between 34 • S and 40 • S, and there were more rapidly latitudinal fluctuations between 40 • S and 55 • S. PM 2.5 concentrations also had an increasing trend (1.02 µg/m 3 ). The global distribution of PM 2.5 was shown in Figure 1c. The air quality in America, Europe, and eastern Russia improved greatly. PM 2.5 concentrations in Africa as a whole showed an increasing trend, and most regions had increased by 10 µg/m 3 . PM 2.5 in Asia, such as Saudi Arabia, India, and China, increased noticeably, indicating serious air pollution. Slight fluctuations were found in Australia.
In the Northern Hemisphere, PM 2.5 concentrations decreased at mid-latitudes, but increased markedly at low latitudes. In the Southern Hemisphere near the tropics, the concentrations of PM 2.5 increased slightly.     Table 2.

Correlation Analysis
The Pearson correlation coefficients were calculated based on each pixel, as shown in Figure 2 (left). The results indicated a significant positive correlation between the LST and SINDVI in most areas of the Northern Hemisphere. In contrast, a highly significant negative correlation appeared in most areas of the Southern Hemisphere, especially Brazilian plateau, Pampas steppe, southeastern Africa, and Australia.

Dominant Factor Analysis
The partial correlation coefficients (i.e., r LST SINDVI (PM2.5) and r LST PM2.5 (SINDVI)) were calculated based on pixels to analyze the main factors causing surface temperature changes (Figure 3). If r LST SINDVI (PM2.5) was larger than r LST PM2.5 (SINDVI), the SINDVI was considered to be the main factor. The opposite result indicated PM2.5 had a larger effect on LST. Latitudinal patterns of Pearson correlation were shown in Figure 2 (right). In the Northern Hemisphere, a positive correlation between LST and SINDVI was found at high latitudes, but a negative correlation was found in other latitudinal zones. In the Southern Hemisphere, the correlation between LST and SINDVI increased then decreased along with a rising latitude. Strong negative correlations between LST and SINDVI occurred at 15-35 • S (i.e., mid-to-low latitudes in the Southern Hemisphere), and the negative correlation gradually weakened from latitude of 40 • S.
Except for low latitudes in the Northern Hemisphere, a positive correlation between LST and PM 2.5 appeared around the world. Compared with regions in north of the equator, those in south of the equator had stronger positive correlations. Positive correlations between LST and PM 2.5 existed in central Eurasia, South America, and eastern Australia, where PM 2.5 fluctuated slightly. The negative correlations between LST and PM 2.5 were mainly concentrated in Europe, West Asia, and south of Sahara Desert. In these areas, PM 2.5 experienced dramatic changes, increasing in West Asia especially in Saudi Arabia and decreasing in Europe such as Poland, Germany, and Romania. Therefore, the effect of PM 2.5 on LST was complex, and it was related to the changes of PM 2.5 concentrations. A small amplitude increasing of PM 2.5 concentration increased LST due to the heat retention effect of PM 2.5 on global surfaces, with a result of a positive correlation. However, a marked increase in PM 2.5 concentration prevented solar radiation from reaching the surface, reducing surface energy and causing a decrease in LST.

Dominant Factor Analysis
The partial correlation coefficients (i.e., r LST SINDVI (PM2.5) and r LST PM2.5 (SINDVI) ) were calculated based on pixels to analyze the main factors causing surface temperature changes ( Figure 3). If r LST SINDVI (PM2.5) was larger than r LST PM2.5 (SINDVI) , the SINDVI was considered to be the main factor. The opposite result indicated PM 2.5 had a larger effect on LST.

Significant Decreases in LST on Regions of Interest (ROIs)
Four regions of interest (ROIs) where LST decreased sharply were found in Saudi Arabia (ROIa), India (ROI-b), Northeast China (ROI-c), and Australia (ROI-d). The spatial positions of ROIs were shown in Figure 1a in yellow and the locations were shown in Table 2. 3.4.1. Saudi Arabia In most parts of the world, SINDVI was found to be the dominant factor (Figure 3). SINDVI increased markedly in the mid-latitude of the Northern Hemisphere, especially over the East European Plain, Labrador Plateau, and India. It decreased obviously in east of the Brazilian Plateau and southeast of Africa. The decrease in vegetation cover caused by human activities led to the increase in LST.

Significant Decreases in LST on Regions of Interest (ROIs)
Four regions of interest (ROIs) where LST decreased sharply were found in Saudi Arabia (ROI-a), India (ROI-b), Northeast China (ROI-c), and Australia (ROI-d). The spatial positions of ROIs were shown in Figure 1a in yellow and the locations were shown in Table 2. 3.4.1. Saudi Arabia LST in Saudi Arabia (ROI-a) has decreased by 2.09 K over the past 16 years (Figure 4a). SINDVI in this region also decreased with a slope of −0.03, but PM 2.5 increased dramatically, with a slope of 1.512. The interannual variations in ROI-a were shown in Figure 4a. The correlation analysis showed a significantly positive correlation between LST and SINDVI (r = 0.661, p < 0.01) and a significantly negative correlation between LST and PM 2.5 (r = −0.746, p < 0.01), which means that both SINDVI and PM 2.5 contributed to the changes in LST. However, the partial correlation analyses indicated lower confidence levels (p > 0.05, Table 3). Therefore, the effects of vegetation cover cannot be ignored in this area.

India
LST experienced a marked decrease as SINDVI and PM 2.5 increased in India (Figure 1). The interannual variations in ROI-b were shown in Figure 4b. The variation in vegetation cover played the leading roles in affecting LST changes during the period from 2001 to 2016 (Figure 3). The correlation analysis showed that LSTs were significantly correlated to both SINDVI and PM 2.5 , with Pearson correlation coefficients of −0.877 (p < 0.01) and −0.630 (p < 0.01), respectively. However, the partial correlation analysis indicated that the correlation between LST and SINDVI (r LST SINDVI (PM2.5) = −0.789, p < 0.01) was more significant than that between LST and PM 2.5 (r LST PM2.5 (SINDVI) = −0.106).

China
LST showed an increasing trend in West China, but a reverse trend in East China, particularly in Northeast China (Figure 1a). The interannual variations in ROI-c were shown in Figure 4c. Interestingly, the partial correlation analyses indicated a significantly positive correlation between LST and SINDVI and a significantly negative correlation between LST and PM 2.5 , while the simple correlation coefficients were smaller than the partial correlation coefficients (Table 3). Therefore, the interaction between SINDVI and PM 2.5 weakened their respective effects on LST in the Northeast China.

Australia
LST in Australia showed a decreasing trend, with a noticeable increase in SINDVI and a slight decrease in PM 2.5 (Figure 1). The interannual variations in ROI-d were shown in Figure 4d. In 2010, the value of LST and PM 2.5 were low, but SINDVI reached its peak (Figure 4d). The simple correlation analysis showed significant relationship between LST with SINDVI and PM 2.5 , but the partial correlation analysis indicated that SINDVI had a larger impact on LST than PM 2.5 (Table 3).

China
LST showed an increasing trend in West China, but a reverse trend in East China, particularly in Northeast China (Figure 1a). The interannual variations in ROI-c were shown in Figure 4c. Interestingly, the partial correlation analyses indicated a significantly positive correlation between  LST is very sensitive to the global climate changes. It had a prominent increasing trend in the Arctic, but a decreasing trend over the Antarctic. This is consistent with the conclusion of Overland et al. [64] that the North and South poles experienced opposite climate changes. Authors identified anthropogenic influence as the primary cause for sea ice reduction in the Arctic [65]. Under the background of global warming, researchers have paid more attention to the warming in the Arctic due to the importance of global sea level rise, especially regarding the Greenland ice sheet [66]. LST could also be used to monitor glacial melt zones based on remote sensing technology. LSTs increased in the Arctic, and the surface temperature of North Greenland had risen by approximately 1 K over the past 16 years (Figure 1a). The main indicators of Arctic temperature increase include reducing sea ice and permafrost. In addition, LST could indicate the burn severity of forest fires [67]. The increase in LST (0.073-0.325 K) has deteriorated forest fires in Siberian boreal forests. In addition to direct warming by fires, lower albedo with darker surfaces and lower evapotranspiration due to biomass burning can increase energy absorption by land surface and decrease the evaporative cooling, with the consequence of higher LST [68].
Although both the Arctic and Antarctic were subjected to similar solar radiation and increasing greenhouse gas concentrations, the two regions had undergone different variations [64]. Studies showed that decreasing LSTs in the Antarctic may be related to the ozone hole [69]. Turner et al. [70] suggested that the ozone hole would shield global warming in the Antarctic. Ozone depletion in the stratosphere had intensified polar vortices, and winds over the center of the Antarctic continent had strengthened the change of weather patterns over the mainland. Strong winds have effectively separated the Antarctic from global warming, with a result of a decrease in LST and an increase in snowfall [71].
SINDVI value was used to denote interannual variations of vegetation productivity [58]. It integrated NDVI over the growth season and weaken the details of NDVI variation. Vegetation cover measured by SINDVI showed a greening trend at the global scale from 2001 to 2016 (Figure 1b). Specifically, vegetation greening trends increased in midlatitude areas of the Northern Hemisphere. Other studies demonstrated that snow cover decreased dramatically in those regions [57]. As a result of global warming, the melting of snow and ice produced a large amount of water that could be used for irrigating plant. In addition, longer growing season caused by more heat can also led to the increases in SINDVI. A dramatic browning trend was found in central Asia (Figure 1b). In this area, a warming climate in the summer worsened the aridity [72]. The decreases in precipitation and aridity were the main reasons for the decreases in SINDVI [72].
The sources of air pollutants were diversified, such as biomass burning, fossil fuel combustion, and industrialization [73,74]. PM 2.5 in the atmosphere could be emitted directly by combustion or mechanical processes and it can also be generated by chemical reactions of inorganic or organic molecules in the atmosphere [73]. The time series analyses of the interannual mean PM 2.5 concentration demonstrated an increasing trend at a global scale (Figure 1c and Table 1). Dramatic decreases in PM 2.5 were observed in the Eastern USA and Europe because of effective pollution control measures, but an opposite trend was found in Sub-Saharan Africa due to biomass burning and higher concentration of mineral dust [60]. New research showed that the emissions of sulfur dioxide had increased by 50% in India from 2007 to 2017 [75]. The sharp increase in PM 2.5 in India was influenced by both anthropogenic and natural factors such as chemical processes [73,75] and mineral dust [60]. The regions where PM 2.5 was the dominate factor for the change in LST were mainly distributed in central Russia and Kazakhstan where the concentration of PM 2.5 increased because of the biomass burning [76], and Saudi Arabia where a lot of mineral dust contributed to the increase in PM 2.5 .

Influence of SINDVI and PM 2.5 on LST
This study analyzed the effects of the vegetation cover measured by SINDVI and air pollution measured by PM 2.5 on the variation in LST. Previous studies showed the lower LSTs with higher vegetation cover [18,77]. Similarly, the current study found a significant negative relationship between SINDVI and LST in the Southern Hemisphere and low latitude areas of the Northern Hemisphere. On the contrary, a positive relationship between LST and SINDVI was found in high latitude areas. Different from the low latitude of Northern Hemisphere, high latitude areas are characterized by lesser solar radiation and lower temperature. Plant photosynthesis is limited particularly by lower temperature and increase in temperature including LST will promote plant growth and NDVI values in high latitude. Therefore, a positive relationship between LST and SINDVI was found in high latitude. A study found that although most regions experienced overall greenness, the increase rate decreased and browning trend expanded in the Northern Hemisphere especially in northern mid-low latitudes [26]. It could be the reason for that the negative relationship between SINDVI and LST in the Northern Hemisphere was weaker than that in the Southern Hemisphere ( Figure 3). Furthermore, compared with past time, the increasing temperature caused a longer plant growth season. It partially caused the increase in SINDVI [78]. Therefore, the browning trend was hidden in the greenness in the Northern Hemisphere. In addition, drought can cause higher LST, but the legacies of drought were different among forest, shrubs, and grass in the North Hemisphere [79].
SUHI and the air pollution have become two major problems along with urbanization without proper environmental protection [43]. The positive correlation between LST and PM 2.5 was found in most areas from 2001 to 2016 (Figure 2b). However, negative correlation was also identified in Saudi Arabia, northeast China and south of Sahara Desert. It might be caused by a different rate of the variation in PM 2.5 concentrations. Cooling effect of aerosol was found in the urban areas [44]. Radiative effects of aerosols, reducing solar shortwave radiation and increasing atmospheric longwave radiation, have a marked cooling effect on air temperature, especially downwind of the city cluster [42]. The cooling effects of aerosol can explain the negative relation between LST and PM 2.5 in Saudi Arabia, northeast China and south of Sahara Desert where PM 2.5 increased markedly. The increase of PM 2.5 can reduce solar radiation reaching land by scattering more radiation back to space, but it could also reduce the surface longwave flux and latent heat flux which may increase LST [80,81].
In addition, Hajiloo et al. [37] used a geographically weighted regression model to assess the relationship between environmental parameters and PM 2.5 concentrations. They found a negative correlation between LST with PM 2.5 in western Tehran and a positive correlation in eastern Tehran which resulted from asymmetry distribution of vegetation. However, different trends of SINDVI did not cause opposite correlation in ROI-a and ROI-b (Figure 4a,b). Therefore, other factors may play a more important part in the negative correlation between LST and PM 2.5 in the two areas.

Limitations and Future Research
Same as many studies, there are limitations in the current research. Due to the spatial heterogeneity, LST can be affected by many factors, such as elevation, albedo, meteorological parameters, human activities, and others. [9,18,57,82]. However, because of the limitation of high spatial resolution remote sensing images and statistical data at a global scale, other factors rather than SINDVI and PM 2. 5 have not yet been considered in this study. In addition, the enhanced vegetation index (EVI) was more sensitive than NDVI in some regions, such as tropical rainforests with high biomass [83]. The saturation of NDVI might underestimate the correlation between LST and vegetation cover. In future studies, other vegetation indexes can be used to analyze vegetation changes in the regions with saturation of NDVI. In addition, it is important to develop models including more climatic, topographical, and socio-economic factors such as precipitation, wind speed, elevation, SUHI, population, and others for better estimate the changes of LST.

Conclusions
In this study, the spatiotemporal trends of global land surface temperature, vegetation cover, and air quality were analyzed during the period of 2001-2016. The research evaluated the influence of vegetation cover and PM 2.5 on land surface temperature. A total of 62.4% of the world's land surface temperature exhibited a increasing trend from 2001 to 2016. A total of 58.2% of regions had an increasing trend of SINDVI, indicating that vegetation coverage increased from a global aspect. The concentrations of PM 2.5 showed an increasing pattern mainly in developing countries, for example China and India, and a decreasing pattern mostly in developed countries, such as the USA. A significant negative correlation between vegetation and land surface temperature indicated vegetation was the main driver for land surface temperature changes in most areas on the earth. The effect of PM 2.5 on surface temperature was not uniform. A significant increase in PM 2.5 concentrations caused a decrease in surface temperature, but a slight increase in PM 2.5 concentrations increased surface temperature. In future studies, models including more parameters should be developed to more accurately analyze the variation in land surface temperature.