Pollution Trends in China from 2000 to 2017: A Multi-Sensor View from Space

: Satellite sensors can provide unique views of global pollution information from space. In particular, a series of aerosol and trace gas monitoring instruments have been operating for more than a decade, providing the opportunity to analyze temporal trends of major pollutants on a large scale. In this study, we integrate aerosol products from MODIS (MODIS Resolution Imaging Spectroradiometer, all abbreviations and their deﬁnitions are listed alphabetically in Abbreviations) and MISR (Multi-angle Imaging Spectroradiometer), the AAI (Absorbing Aerosol Index) product from OMI (Ozone Monitoring Instrument), column SO 2 and NO 2 concentrations from OMI, and tropospheric column ozone concentration from OMI / MLS (Microwave Limb Sounder) to study temporal changes in major pollutants over China. MODIS and MISR consistently revealed that column AOD (Aerosol Optical Depth) increased from 2000, peaked around 2007, and started to decline afterward, except for northwest and northeast China, where a continuous upward trend was found. Extensive negative trends in both SO 2 and NO 2 have also been found over major pollution source regions since ~2005. On the other hand, the OMI AAI exhibited signiﬁcant increases over north China, especially the northeast and northwest regions. These places also have a decreased Angstrom exponent as revealed by MISR, indicating an increased fraction of large particles. In general, summer had the largest AOD, SO 2 , and NO 2 trends, whereas AAI trends were strongest for autumn and winter. A multi-regression analysis showed that much of the AOD variance over major pollution source regions could be explained by SO 2 , NO 2 , and AAI combined, and that the SO 2 and NO 2 reduction was likely to be responsible for the negative AOD trends, while the AOD increase over NE and NW China may be associated with an increase of coarse particles revealed by increased AAI and decreased AE. In contrast to aerosols, tropospheric ozone exhibited a steady increase from 2005 throughout China. This indicates that although the recent emission control e ﬀ ectively reduced aerosol pollutants, ozone remains a challenging issue and may dominate future air pollution.


Introduction
In recent years, China has been facing a challenging air pollution problem in the background of its rapid industrialization and economic development. Especially in the past two decades, the frequent occurrence of haze and smog days has brought about considerable damage to both the environment and human health [1] as well as perturbations to the climate [2]. As a result, pollution in China has become the focus of many studies in China and worldwide. The major emission source of anthropogenic pollutants in China is fossil fuel burning, especially coal burning [3,4]. Studies have shown that this source has contributed to over 80% of SO 2 pollution in China [5]. In urban areas, traffic is also an important source that greatly contributes to NO x , ozone, and particulate matter pollution [3,6,7]. In the winter, north China typically suffers from increased black carbon emissions from heating [8,9]. Seasonal biomass burning [10] and dust storms [11] also impair air quality over north and southeast The MODIS instrument is a single-view imager with a swath width of 2330 km and near global coverage every two days. This high sampling frequency captures most of the aerosol variability and microphysics properties. MODIS was launched on EOS Terra in 1999 and on EOS Aqua in 2002. We used Collection 6.1 Level 2 AOD products [26,27] from both of these platforms. The variable used was the quality controlled dark target and deep blue combined AOD at 550 nm over land, and the time period spanned from 2000 to 2017. The land AOD product uncertainty was estimated to be 0.05 ± 20%AERONET (Aerosol Robotic Network) AOD. Only QC = 3 data were selected to ensure data quality. We gridded the pixel level data into daily means at 0.5 • × 0.5 • resolution, and only grids with more than half the available pixel level data were considered valid. To calculate the monthly means, we applied two schemes, namely threshold equal day averaging and pixel weighted averaging, as Levy et al. [28] indicated that substantial differences in the monthly means could be caused by different averaging schemes. The two schemes were found to have a negligible impact on the trend (figure not Remote Sens. 2020, 12, 208 3 of 24 shown). We thus only show the results using 5-day threshold equal day averaging (i.e., each month must have at least five valid daily means).
MISR is a multi-angle sensor with nine pushbroom cameras on the EOS Terra platform. The zonal overlap of the common swath of all nine cameras is at least 360 km in order to provide multi-angle coverage in nine days at the equator, and two days at the poles [29]. Compared to MODIS, the multi-angle view of MISR performs better over bright surfaces, while its lower sampling may not fully resolve small-scale variability. In addition to AOD, the AE is also a useful parameter indicating the composition of fine/coarse particles. Therefore, we used version 31 Level 3 gridded monthly mean AOD and AE products [30][31][32] with a 0.5 • resolution from 2000 to 2017. The reported MISR AOD uncertainty is 0.05 or 20%AOD [33,34]. Note that MODIS AE over land is considered unsuitable for quantitative research [13] and we thus only used those data from MISR.

OMI Absorbing Aerosol Index
The OMI launched onboard the EOS Aura in 2004 [35] provides an absorbing aerosol index (AAI) product, which constituted a qualitative indicator of the presence of UV absorbing aerosols in the atmosphere such as biomass burning and dust [36]. The AAI calculations were based on the difference in surface reflectivity derived from the 331.2 and 360 nm backscatter radiance, exhibiting an uncertainty of ±30% [36]. The AAI dataset used in this study was the version 3 Level 2 daily product with 0.25 • resolution from the near UV algorithm [37,38] during the 2005-2017 period. As clouds and scattering aerosols can also produce small positive AAI, we only retained AAI values greater than 0.7 to represent absorbing aerosols as suggested by Prospero et al. (2002) [39]. Then, the monthly means were calculated with a 5-day minimum threshold. Although more quantitative products such as the absorption aerosol optical depth and single scattering albedo are also provided, these data tend to have greater uncertainty arising from aerosol model assumptions in the retrieval algorithm. Therefore, we used the AAI data here, which is the ratio of the directly observed radiance at two channels.

OMI SO 2 and NO 2
The measurements of SO 2 and NO 2 are also explicit objectives of OMI. We used the operational OMI planetary boundary layer SO 2 product processed with the improved band residual difference algorithm [40][41][42]. The algorithm uses differential residuals at the three wavelength pairs with the largest differential SO 2 cross sections in the OMI UV2 spectral region (310-380 nm) to maximize measurement sensitivity to anthropogenic emissions in the PBL (Planetary Boundary Layer). Validation against field measurements indicates accuracy of~1.5 DU (Dobson units) over East China [43].
The tropospheric column NO 2 data used are from the NASA standard product version 3 [41]. It uses the new algorithm described in [43]. The average error reported in the data file was 2.2 × 10 15 molecules/cm 2 . For consistency with SO 2 , the data were converted from molecules/cm 2 to Dobson unit (1 DU = 2.69 × 10 16 molecules/cm 2 ).
The original product used was version 3, Level 3 daily mean SO 2 and NO 2 products with 0.25 • resolution [38]. The data were further gridded into monthly means at 0.5 • × 0.5 • resolution, for the purpose of multi-regression analysis. The study period was also from 2005 to 2017.

OMI/MLS Tropospheric Ozone
As OMI only retrieves the total column ozone, a joint algorithm was developed to retrieve the tropospheric ozone with the MLS also onboard Aura [44]. The OMI ozone data were provided on a 1 • × 1.25 • grid for each day. There are typically 3494 MLS Level 2 Version 1.5 along track profile per day over the globe with horizontal resolution of 150 km and vertical resolution of 3 km in the upper troposphere and lower stratosphere. OMI is matched with each MLS profile by linearly interpolating the nearest four 1 • × 1.25 • grids. Then, the corresponding MLS ozone profile was integrated from the tropopause to the top of the atmosphere and extracted from the matched OMI total column ozone. The final tropospheric ozone data was gridded with 1 • × 1.25 • grid for each month [45]. The OMI/MLS tropospheric ozone data compared well with ozonesondes and chemical transport models [46,47] in terms of large scale features including the mean value, seasonal cycle, and zonal variations, although over high latitudes and the Sahara desert, a larger different exceeding 5 DU can be observed.

Trend Estimation
The linear trends are estimated using the monthly mean data with the Sen's slope method [48]. This is a non-parametric estimation and is less sensitive to outliers compared to linear regression. The Sen's slop is defined as the median of the slopes of all data pairs: where b is the Sen's slope and Xi and Xj are the ith and jth data in the time series, respectively.
To test the significance of the trend, we used the Mann-Kendall statistical test [49,50]. This is also a nonparametric test to identify the existence of monotonic trends in a time series. The MK test statistic is calculated as: which follows a normal distribution with zero mean and variance.
Moreover, studies have found that the existence of serial correlation in the time series will increase the probability that the MK test detects a significant trend [51][52][53]. We thus pre-whitened the time series to remove serial correlation following the procedure recommended by [53]. Specifically, the linear trend is first removed from the time series using the Sen's slope b: Then the AR(1) process is estimated and removed: where r 1 is the lag−1 autocorrelation, and finally the linear trend is added back:

Multiple Regression
To investigate the relationship between AOD and the components (AAI, SO 2 , and NO 2 as precursors for sulfate and nitrate aerosols respectively), we performed a multiple regression analysis as: Remote Sens. 2020, 12, 208

of 24
The time series for AOD, AAI, SO 2 , and NO 2 are first de-seasonalized and de-trended, by removing the multi-year averaged seasonal cycle and linear trends, respectively, and then normalized by their respective standard deviation. In this way, the magnitude of each regression coefficient can reflect the relative contribution of the corresponding variable to the total variance of AOD. The significance of the model as well as each coefficient is tested using the F test. In addition, because significant correlation may exist between some species such as SO 2 and NO 2 , which will impact on the accuracy of the multi-regression coefficients [54,55], we also performed a single parameter regression of AOD against each species to verify the credibility of the coefficients.

Results
Before discussing the trends, we first briefly examine the spatial distribution of the parameters used in this study ( Figure 1). It was seen that the highest AOD values (up to 1) mainly concentrated on East China and the Sichuan Basin, corresponding to the most densely populated areas. NW China's Taklimakan Desert also showed some high AOD signals, resulting from wind-blown dust. East China also had the highest SO 2 and NO 2 concentrations, especially the NCP (North China Plain) region. AE exhibited higher values over East China than the west, and the lowest values were found in the Taklimakan Desert. This is also consistent with the general aerosol size distribution as East China is dominated by fine anthropogenic aerosols whereas dust is typically found over the desert region in the west. The spatial distribution of AAI also showed a hotspot in the Taklimakan Desert, with some moderately high values found over the NCP and NE China. The former is the result of dust aerosols, while the latter should be mainly attributed to carbonaceous aerosols considering the higher AE values there. The spatial distribution of tropospheric O 3 was more spatially uniform than the aerosol and the other two gas components, and largely follows the topography (i.e., higher concentrations over low altitudes and lowest concentration over the Tibetan Plateau). In general, the parameters examined reflect the spatial distribution of major pollutants in China.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 25 The time series for AOD, AAI, SO2, and NO2 are first de-seasonalized and de-trended, by removing the multi-year averaged seasonal cycle and linear trends, respectively, and then normalized by their respective standard deviation. In this way, the magnitude of each regression coefficient can reflect the relative contribution of the corresponding variable to the total variance of AOD. The significance of the model as well as each coefficient is tested using the F test. In addition, because significant correlation may exist between some species such as SO2 and NO2, which will impact on the accuracy of the multi-regression coefficients [54,55], we also performed a single parameter regression of AOD against each species to verify the credibility of the coefficients.

Results
Before discussing the trends, we first briefly examine the spatial distribution of the parameters used in this study ( Figure 1). It was seen that the highest AOD values (up to 1) mainly concentrated on East China and the Sichuan Basin, corresponding to the most densely populated areas. NW China's Taklimakan Desert also showed some high AOD signals, resulting from wind-blown dust. East China also had the highest SO2 and NO2 concentrations, especially the NCP (North China Plain) region. AE exhibited higher values over East China than the west, and the lowest values were found in the Taklimakan Desert. This is also consistent with the general aerosol size distribution as East China is dominated by fine anthropogenic aerosols whereas dust is typically found over the desert region in the west. The spatial distribution of AAI also showed a hotspot in the Taklimakan Desert, with some moderately high values found over the NCP and NE China. The former is the result of dust aerosols, while the latter should be mainly attributed to carbonaceous aerosols considering the higher AE values there. The spatial distribution of tropospheric O3 was more spatially uniform than the aerosol and the other two gas components, and largely follows the topography (i.e., higher concentrations over low altitudes and lowest concentration over the Tibetan Plateau). In general, the parameters examined reflect the spatial distribution of major pollutants in China.

AOD Trends
The linear trends for MODIS and MISR AOD are shown in Figure 2

AOD Trends
The linear trends for MODIS and MISR AOD are shown in Figure 2. During the whole period of 2000 to 2017, the MODIS Terra data indicated weak increasing trends over the NCP and NE China and weak decreasing trends over Inner Mongolia to the northwest of the NCP (Figure 2a), MODIS Aqua exhibited significant negative trends over South China and insignificant positive trends over NE China (Figure 2d), while MISR showed weak significant negative trends in parts of South China, but largely insignificant trends elsewhere (Figure 2g). Through closer examination, we found that there appeared to be a sharp trend change in 2008. Before 2008, both MODIS and MISR data showed strong and significant positive changes as large as 0.3/decade over the entire East China (middle column). After 2008, in contrast, prevailing negative trends with comparable magnitudes were found over China (right column). The only exception is in NE and NW China, where positive trends were noticed. This trend shift around 2008 was also previously noticed by [56] using MODIS data.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 25 Aqua exhibited significant negative trends over South China and insignificant positive trends over NE China (Figure 2d), while MISR showed weak significant negative trends in parts of South China, but largely insignificant trends elsewhere (Figure 2g). Through closer examination, we found that there appeared to be a sharp trend change in 2008. Before 2008, both MODIS and MISR data showed strong and significant positive changes as large as 0.3/decade over the entire East China (middle column). After 2008, in contrast, prevailing negative trends with comparable magnitudes were found over China (right column). The only exception is in NE and NW China, where positive trends were noticed. This trend shift around 2008 was also previously noticed by [56] using MODIS data. To corroborate satellite retrieval products, we also examined the trends at two ground sites within the AERONET, which is a surface network of high quality sunphotometer measurements of aerosol properties, whose data are often used to validate satellite retrievals [57]. The two sites were Beijing and XiangHe, which are both located in the NCP and are the only two sites with sufficiently long data records for long term trend estimation. From Figure 3, it is clearly seen that both these sites exhibited trend shifts from positive to negative around 2008, which is consistent with the satellite data. To corroborate satellite retrieval products, we also examined the trends at two ground sites within the AERONET, which is a surface network of high quality sunphotometer measurements of aerosol properties, whose data are often used to validate satellite retrievals [57]. The two sites were Beijing and XiangHe, which are both located in the NCP and are the only two sites with sufficiently long data records for long term trend estimation. From Figure 3, it is clearly seen that both these sites exhibited trend shifts from positive to negative around 2008, which is consistent with the satellite data.
The seasonal AOD trends from MODIS showed a very similar pattern to the annual trends (Figures 4 and 5) and little differences were found between the results of the four seasons. To save space, only the MODIS Terra results are presented here because the Aqua results were quite similar. In general, before 2008, the winter (DJF, December, January, February) showed the strongest increase, especially over the NCP and PRD, where the trends exceeded 0.3/decade, whereas the trend during the other seasons was generally between 0.1-0.2/decade. After 2008, the summer (JJA, June, July, August) and fall (SON, September, October, November) trends appeared slightly stronger, especially for the southern part. Note that we did not use MISR for seasonal analysis, mainly because we observed more variability and missing data in the MISR time series, probably related to its lower sampling frequency than MODIS. This sampling difference may also explain many of the trend differences between MODIS and MISR in Figure 2. The seasonal AOD trends from MODIS showed a very similar pattern to the annual trends (Figures 4 and 5) and little differences were found between the results of the four seasons. To save space, only the MODIS Terra results are presented here because the Aqua results were quite similar. In general, before 2008, the winter (DJF, December, January, February) showed the strongest increase, especially over the NCP and PRD, where the trends exceeded 0.3/decade, whereas the trend during the other seasons was generally between 0.1-0.2/decade. After 2008, the summer (JJA, June, July, August) and fall (SON, September, October, November) trends appeared slightly stronger, especially for the southern part. Note that we did not use MISR for seasonal analysis, mainly because we observed more variability and missing data in the MISR time series, probably related to its lower sampling frequency than MODIS. This sampling difference may also explain many of the trend differences between MODIS and MISR in Figure 2.

AAI and AE Trends
The AAI is defined as the difference between the ratio of two satellite observed radiances at the UV channel and those calculated for a pure Rayleigh atmosphere. It is a measure of the wavelength-dependent reduction of Rayleigh scattered radiance by aerosol absorption. By theory, it is related to the amount and height of UV absorbing aerosols [23] including dust, black, and organic carbon. Therefore, it offers insights into aerosol composition. Figure 6a shows the annual trends of OMI AAI during the 2005 to 2017 period (the trends from 2008 to 2017 are similar). Significant increases were noted over most of North China. The strongest positive trends were found in the Taklimakan Desert in the Xinjiang Province of NW China. The trends there reached~0.27/decade, while over the other regions with positive trends, the trend amplitude was mostly below 0.2/decade. For the southern part where significant negative AOD trends were observed during this period, the AAI trend was weak and insignificant.
Since the AAI increases can be caused by changes in the loading of either dust or carbonaceous aerosols or both, we also examined the trend for MISR AE (Figure 6b). The AE is related to particle size and a lower AE indicates larger sizes in general. Comparing Figures 6a and 6b, it can be clearly seen that the regions with positive AAI changes corresponded well to those with negative AE trends. This suggests an increase in the mean particle size, which is mostly likely due to increased dust fraction. Considering that NW China is also a major dust source region, increased dust activity can well explain the observed AAI and AE trends. The trends in NE China might be associated with increased dust transport from the west.
Seasonally, the strongest AAI increase happened in the fall and winter ( Figure 7). This is somewhat inconsistent with the dust outburst season, which is usually in the spring. However, precious studies have also indicated increased dust emission in the late fall-winter season [58,59]. We also noticed the strongest negative trends in MISR AE for SON (figure not shown), while for winter, the MISR data availability was too low to estimate significant trends. size and a lower AE indicates larger sizes in general. Comparing Figure 6a and Figure 6b, it can be clearly seen that the regions with positive AAI changes corresponded well to those with negative AE trends. This suggests an increase in the mean particle size, which is mostly likely due to increased dust fraction. Considering that NW China is also a major dust source region, increased dust activity can well explain the observed AAI and AE trends. The trends in NE China might be associated with increased dust transport from the west. Seasonally, the strongest AAI increase happened in the fall and winter (Figure 7). This is somewhat inconsistent with the dust outburst season, which is usually in the spring. However, precious studies have also indicated increased dust emission in the late fall-winter season [58,59]. We also noticed the strongest negative trends in MISR AE for SON (figure not shown), while for winter, the MISR data availability was too low to estimate significant trends.
Relating the AAI trends to those of AOD, it was noticed that the regions with an AAI increase largely corresponded to regions with positive AOD trends during the 2008-2017 period (i.e., NE and NW China, as seen in Figure 2c). The seasonality was also in phase with the most notable changes in the fall and winter, although the AOD trends appeared weaker. This makes sense since dust only accounts for a small fraction of total aerosols in these two seasons. Moreover, since AAI is not linearly related to dust concentration, its change cannot reflect the absolute change in the dust amount.  It is worth noting that AAI is a qualitative indicator of absorbing aerosols and its trends should be interpreted with caution. According to previous studies [34,60,61], apart from the column loading of absorbing aerosols, AAI is also sensitive to absorbing aerosol layer height, surface albedo, clouds, and ozone concentration. We then discuss the possible impact of these four factors on the observed trends in detail. First, the aerosol layer height is typically highly correlated with PBLH (Planetary Boundary Layer), and previous studies have indicated that AAI is positively related to the PBLH. However, Guo et al. (2019) [62] revealed a downward PBLH trend over China since 2004, which means that the change of BLH would cause a decreasing AAI trend, opposite to what we have shown. Next, the increase of surface albedo has two competing effects on AAI: a higher albedo may enhance aerosol absorption and thus increase AAI, while also increasing the total amount of reflected radiation at the top of the atmosphere, which may reduce the AAI. De Graaf et al. (2005) [60] indicated that the overall effect will be dominated by the first mechanism (i.e., increasing surface albedo may increase AAI). Although a long-term surface albedo dataset is not readily available, a recent study by Relating the AAI trends to those of AOD, it was noticed that the regions with an AAI increase largely corresponded to regions with positive AOD trends during the 2008-2017 period (i.e., NE and NW China, as seen in Figure 2c). The seasonality was also in phase with the most notable changes in the fall and winter, although the AOD trends appeared weaker. This makes sense since dust only accounts for a small fraction of total aerosols in these two seasons. Moreover, since AAI is not linearly related to dust concentration, its change cannot reflect the absolute change in the dust amount.
It is worth noting that AAI is a qualitative indicator of absorbing aerosols and its trends should be interpreted with caution. According to previous studies [34,60,61], apart from the column loading of absorbing aerosols, AAI is also sensitive to absorbing aerosol layer height, surface albedo, clouds, and ozone concentration. We then discuss the possible impact of these four factors on the observed trends in detail. First, the aerosol layer height is typically highly correlated with PBLH (Planetary Boundary Layer), and previous studies have indicated that AAI is positively related to the PBLH. However, Guo et al. (2019) [62] revealed a downward PBLH trend over China since 2004, which means that the change of BLH would cause a decreasing AAI trend, opposite to what we have shown. Next, the increase of surface albedo has two competing effects on AAI: a higher albedo may enhance aerosol absorption and thus increase AAI, while also increasing the total amount of reflected radiation at the top of the atmosphere, which may reduce the AAI. De Graaf et al. (2005) [60] indicated that the overall effect will be dominated by the first mechanism (i.e., increasing surface albedo may increase AAI). Although a long-term surface albedo dataset is not readily available, a recent study by Chen et al. (2019) [63] found that the leaf area of vegetation has significantly increased over China since 2000 by using satellite measurements. This suggests that surface albedo is likely to be decreased during this period as vegetation has a lower albedo than bare lands. This decreased surface albedo would lead to a decreased AAI trend, which is again in contrast to the observation. Third, the effect of clouds on AAI is more complicated and depends on the relative height of the cloud and aerosol layer. The accurate analysis of this effect requires long term observations of the collocated cloud height and aerosol height, which is rather difficult to obtain. Nonetheless, Li et al. (2018) [64] found that cloud fraction over China did not exhibit significant trends, therefore, we consider that the change of clouds is unlikely to be a dominant factor for the AAI change. Finally, absorption of ozone in the UV range will also increase the AAI value. The AAI product already increased the correction of ozone absorption. For the study period, although we noticed an increase in tropospheric ozone over China (Section 3.5), ozone absorption was dominated by stratospheric ozone and the change in tropospheric ozone was only a minor contribution to total ozone absorption. Moreover, the spatial distribution of the ozone trends was not consistent with the AAI trends. These facts suggest that ozone change is also unlikely to be responsible for the AAI trends. In sum, factors other than absorbing aerosol loading that affect the AAI have mostly been ruled out. Combined with the negative AE trends, the positive AAI trends are thus mostly likely to be associated with increases in the amount of absorbing aerosols, especially dust.

SO 2 and NO 2 Trends
SO 2 and NO 2 are also major pollutants in China as well as precursors for sulfate and nitrate aerosols. Their trends directly reflect the emission changes, which may help explain the AOD trends. In Figure 8, we plot the annual trends for the column SO 2 and NO 2 concentrations retrieved by OMI for the 2005-2017 period. The spatial pattern of SO 2 trends highly resembles that of AOD (Figure 2c,f,i), with significant decreases over East and South China, in particular, the major pollution source regions of the NCP, Sichuan Basin, and PRD. The distribution of trends for NO 2 appears to be comparatively scattered and also exhibits significant negative trends over East China, but are mainly concentrated on the southern NCP, YRD, and the PRD. The major source of NO 2 is automotive, so the highest emissions are typically found over densely polluted city conglomerates. Nonetheless, the SO 2 and NO 2 trends both indicate significantly decreased emissions of anthropogenic pollution and well agree with the negative trends of AOD for this time period.
regions of the NCP, Sichuan Basin, and PRD. The distribution of trends for NO2 appears to be comparatively scattered and also exhibits significant negative trends over East China, but are mainly concentrated on the southern NCP, YRD, and the PRD. The major source of NO2 is automotive, so the highest emissions are typically found over densely polluted city conglomerates. Nonetheless, the SO2 and NO2 trends both indicate significantly decreased emissions of anthropogenic pollution and well agree with the negative trends of AOD for this time period. The SO2 trends for the four seasons were highly consistent with the annual trends ( Figure 9). Regions with negative trends extended slightly to the northeast in the summer and fall. For seasonal NO2 trends, some differences were observed ( Figure 10). Specifically, while overall, most areas over East China showed negative trends, during the winter and fall, significant increases were found in the northeast. This could also partly contribute to the increased AOD over NE China during these two seasons, in addition to dust. The SO 2 trends for the four seasons were highly consistent with the annual trends ( Figure 9). Regions with negative trends extended slightly to the northeast in the summer and fall. For seasonal NO 2 trends, some differences were observed ( Figure 10). Specifically, while overall, most areas over East China showed negative trends, during the winter and fall, significant increases were found in the northeast. This could also partly contribute to the increased AOD over NE China during these two seasons, in addition to dust.     The decreases of SO 2 and NO 2 are consistent with previous study by [65], who analyzed global trends of SO 2 and NO 2 from OMI, and are clear indications of the effectiveness of emission reduction measures taken in China during the past decade.

Multiple Regression
Multiple regression was performed on AOD against AAI, SO 2 , and NO 2 according to Equation (7) to further investigate the possible compositional changes associated with the AOD trends. The variance of AOD explained by the multiple regression model, defined as the R 2 value between the original AOD time series and the regressed time series, is shown in Figure 11. Note that the AOD variance here refers to the deviation from its multi-year averaged seasonal cycle. We notice that for the majority of places over North and East China, the regression model was significant and the variance explained reached 50% (corresponding to~0.7 correlation) or higher. The highest variance explained was found over NW China's Taklimakan Desert. This is not surprising because this region is dominated by dust aerosols, which is a primary aerosol species and can be well represented by the AAI. For most of the other regions, especially those in East China, the aerosol composition is more complex. Moreover, many aerosol species are secondary such as sulfate, nitrate, and organic aerosols, thus making the relationship between AOD and gas precursors ambiguous and variable. The regression model still performed well over major pollution source regions such as the NCP, PRD, YRD, and the Sichuan Basin. This is an indication of the success of the fitting, and that the decomposition of AOD variability into that of AAI, SO 2 , and NO 2 is reasonable.
Remote Sens. 2020, 12, 208 13 of 24 AAI. For most of the other regions, especially those in East China, the aerosol composition is more complex. Moreover, many aerosol species are secondary such as sulfate, nitrate, and organic aerosols, thus making the relationship between AOD and gas precursors ambiguous and variable. The regression model still performed well over major pollution source regions such as the NCP, PRD, YRD, and the Sichuan Basin. This is an indication of the success of the fitting, and that the decomposition of AOD variability into that of AAI, SO2, and NO2 is reasonable. Figure 11. Variance of AOD explained by the multiple regression model against AAI, SO2, and NO2. The significance of the regression was tested using the F test and places passing the 95% significance level are marked by black dots.
We continued to examine the regression coefficients as they reflect the relative contribution of each component to the total AOD given that all of the time series were normalized. Figure 12 shows the spatial distribution of the coefficients for AAI, SO2, and NO2. To examine the possible effect of multicollinearity on the accuracy of the regression coefficients, we compared the multi-regression coefficients with those obtained from the single parameter regression of AOD onto AAI, SO2, and NO2, respectively, and found that the magnitude and distribution of the coefficients were highly Figure 11. Variance of AOD explained by the multiple regression model against AAI, SO 2 , and NO 2 . The significance of the regression was tested using the F test and places passing the 95% significance level are marked by black dots.
We continued to examine the regression coefficients as they reflect the relative contribution of each component to the total AOD given that all of the time series were normalized. Figure 12 shows the spatial distribution of the coefficients for AAI, SO 2 , and NO 2 . To examine the possible effect of multicollinearity on the accuracy of the regression coefficients, we compared the multi-regression coefficients with those obtained from the single parameter regression of AOD onto AAI, SO 2 , and NO 2 , respectively, and found that the magnitude and distribution of the coefficients were highly consistent (figure not shown). It can be clearly seen that the signals of AAI coefficients (Figure 12a) were mostly concentrated in the north, especially the desert area of the northwest, whereas the highest coefficients for SO 2 (Figure 12b) and NO 2 (Figure 12c) appeared in the east and the north. Again, this was expected, as dust is the major aerosol component in the northwest, thus high signals of AAI were found there. On the other hand, anthropogenic aerosols are the major contributor to AOD over East and South China, thus the AOD variability of the latter is more closely related to sulfate and nitrate precursors of SO 2 and NO 2 , respectively. Note that even for SO 2 and NO 2 , the distributions of the coefficients were not the same. The SO 2 pattern showed two distinct centers over the NCP and the Sichuan Basin, respectively, with negligible signals elsewhere. For NO 2 , large positive coefficients were found over the entire East Asia and even the northeast, in particular, the coastal regions of Southeast and South China where the SO 2 impact is low. Although SO 2 and NO 2 can both be produced by fossil fuel burning, they have different emission sources. SO 2 is primarily released from coal burning while automobiles release the largest fraction of anthropogenic NO 2 . A large fraction of NO 2 is also formed from natural processes such as lightning. The difference in the spatial distribution of SO 2 and NO 2 regression coefficients may reflect the different aerosol compositions in East and South China. coastal regions of Southeast and South China where the SO2 impact is low. Although SO2 and NO2 can both be produced by fossil fuel burning, they have different emission sources. SO2 is primarily released from coal burning while automobiles release the largest fraction of anthropogenic NO2. A large fraction of NO2 is also formed from natural processes such as lightning. The difference in the spatial distribution of SO2 and NO2 regression coefficients may reflect the different aerosol compositions in East and South China. Seasonal multiple regression was also performed as the residence time and dominant pollution species tend to vary with season. The seasonal models also showed satisfactory fitting with the variance explained mostly over 50% ( Figure 13). The spring season (MAM) had the best overall fitting. The missing values in winter (DJF) over the high latitudes and the Tibetan Plateau were due to the high surface albedo caused by the ice and snow surface. Figure 14 displays the seasonal regression coefficients for the three parameters. For the AAI, large positive signals exist over the Taklamakan Desert in all seasons. In the spring, positive signals extended over most of China, indicating wide spread dust events because this is the windiest season of the year. The positive signals can also extend to the NCP and parts of South China, except for summer (JJA). For the summer season, the concentrations of both dust and absorbing carbonaceous aerosols were in general lower than the other three seasons, which was due to decreased dust emission and increased scattering associated with the photochemical production of secondary species and higher aerosol water uptake. The regression coefficients of SO2 and NO2 showed less seasonality, both agreeing with annual results (Figure 12b,c) and concentrating over densely populated regions including the NCP, Southeast, and South China. Seasonal multiple regression was also performed as the residence time and dominant pollution species tend to vary with season. The seasonal models also showed satisfactory fitting with the variance explained mostly over 50% ( Figure 13). The spring season (MAM) had the best overall fitting. The missing values in winter (DJF) over the high latitudes and the Tibetan Plateau were due to the high surface albedo caused by the ice and snow surface. Figure 14 displays the seasonal regression coefficients for the three parameters. For the AAI, large positive signals exist over the Taklamakan Desert in all seasons. In the spring, positive signals extended over most of China, indicating wide spread dust events because this is the windiest season of the year. The positive signals can also extend to the NCP and parts of South China, except for summer (JJA). For the summer season, the concentrations of both dust and absorbing carbonaceous aerosols were in general lower than the other three seasons, which was due to decreased dust emission and increased scattering associated with the photochemical production of secondary species and higher aerosol water uptake. The regression coefficients of SO 2 and NO 2 showed less seasonality, both agreeing with annual results (Figure 12b,c) and concentrating over densely populated regions including the NCP, Southeast, and South China.  Overall, the multiple regression model using AAI, SO2, and NO2 as independent variables well fit the AOD variability and the spatial distribution of the coefficients was reasonable. Through this, combined with the AAI, SO2, and NO2 trends revealed previously, we can conclude that the decrease of AOD after ~2008 over East and South China is likely to be the result of reduced SO2 and NO2 emissions. Increase of dust (reflected by AAI) should be primarily responsible for the weak positive AOD trends over NW and NE China. However, caution should be taken in interpreting the multiple regression results, because (1) satellite retrievals bear certain sources of uncertainties, (2) the AAI, SO2, and NO2 are not directly related to the amount of dust, sulfate, and nitrate aerosols, and (3) like most statistical methods, correlation does not necessarily mean causality, and coincident measurements as well as modeling study is likely needed to confirm the causal relationship.

Tropospheric Ozone Trends
In contrast to aerosols, tropospheric ozone over China exhibits a nation-wide significant increasing trend on the order of 5 DU/decade ( Figure 15). The spatial distribution of the trends is also much more uniform compared to aerosols, which is also consistent with the more uniform spatial Overall, the multiple regression model using AAI, SO 2 , and NO 2 as independent variables well fit the AOD variability and the spatial distribution of the coefficients was reasonable. Through this, combined with the AAI, SO 2 , and NO 2 trends revealed previously, we can conclude that the decrease of AOD after~2008 over East and South China is likely to be the result of reduced SO 2 and NO 2 emissions. Increase of dust (reflected by AAI) should be primarily responsible for the weak positive AOD trends over NW and NE China. However, caution should be taken in interpreting the multiple regression results, because (1) satellite retrievals bear certain sources of uncertainties, (2) the AAI, SO 2 , and NO 2 are not directly related to the amount of dust, sulfate, and nitrate aerosols, and (3) like most statistical methods, correlation does not necessarily mean causality, and coincident measurements as well as modeling study is likely needed to confirm the causal relationship.

Tropospheric Ozone Trends
In contrast to aerosols, tropospheric ozone over China exhibits a nation-wide significant increasing trend on the order of 5 DU/decade ( Figure 15). The spatial distribution of the trends is also much more uniform compared to aerosols, which is also consistent with the more uniform spatial distribution of tropospheric ozone (Figure 1). This should be reasonable considering the longer lifetime of ozone. The southwest region and edge of the Tibetan Plateau showed slightly higher trends, which is likely to be associated with transport from India and Southeast Asia [13][14][15][16]63] as the trends there appear stronger than in China (not shown in Figure 15). Distinct differences were observed in the seasonal trends ( Figure 16). In general, spring trends were the strongest (5 DU/decade on average) while summer had the most significant trends (over 80% of the grids showed significant trends). The fall (SON) and winter (DJF) trends were mostly below 4 DU/decade and only~40% of the grids showed significant trends. Tropospheric ozone pollution is typically heavier in spring and summer due to the high temperature and stronger surface solar radiation [66]. Therefore, it is natural that these two seasons exhibit greater changes than the other two seasons. This is most prominent for North China, whereas for a few regions including the Tibetan Plateau and Southwest China, larger increases were observed in the fall and winter.
The increase of tropospheric ozone in China has also been reported at scattered sites using ground observations [25,67,68]. Li et al. [69] further attributed this positive trend to the reduction of particular matter. Compared to ground observations, satellites have the advantage of revealing large scale patterns, and help to identify sources and transport features. The results shown here indicate that the worsening of ozone pollution is not local, but rather extends all over China.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 25 distribution of tropospheric ozone (Figure 1). This should be reasonable considering the longer lifetime of ozone. The southwest region and edge of the Tibetan Plateau showed slightly higher trends, which is likely to be associated with transport from India and Southeast Asia [13][14][15][16]63] as the trends there appear stronger than in China (not shown in Figure 15). Distinct differences were observed in the seasonal trends ( Figure 16). In general, spring trends were the strongest (5 DU/decade on average) while summer had the most significant trends (over 80% of the grids showed significant trends). The fall (SON) and winter (DJF) trends were mostly below 4 DU/decade and only ~40% of the grids showed significant trends. Tropospheric ozone pollution is typically heavier in spring and summer due to the high temperature and stronger surface solar radiation [66]. Therefore, it is natural that these two seasons exhibit greater changes than the other two seasons. This is most prominent for North China, whereas for a few regions including the Tibetan Plateau and Southwest China, larger increases were observed in the fall and winter.   The increase of tropospheric ozone in China has also been reported at scattered sites using ground observations [25,67,68]. Li et al. [69] further attributed this positive trend to the reduction of particular matter. Compared to ground observations, satellites have the advantage of revealing large scale patterns, and help to identify sources and transport features. The results shown here indicate that the worsening of ozone pollution is not local, but rather extends all over China.

Regional Trends
To better examine the variability of the pollutants as well as regional trend differences, in Figure  17, we plotted the annual time series of all pollutants over six representative regions marked on Figure 1a. The thick black lines in Figure 17 show the four-year moving average smoothed time series. The percentage of change is also marked on Figure 17 with positive changes shown in red and negative in blue.
The changes of the time series in Figure 17 can mostly be inferred from the trend maps shown previously, but some detailed regional features are still observed. For example, although in general, the AOD has started to decrease since 2008, at NCP, the AOD increased until ~2012 before a declining trend started. NO2 showed similar behaviors for the NCP, YRD, and Sichuan Basin, where it first increased until ~2012 and then experienced a sharp decline. Moreover, different to most of the country where overall negative AOD, NO2, and SO2 trends in the last decade were found, over NW China, these parameters were steadily increasing, together with AAI. This may be partly associated with the China western development strategy, which increased the emission of anthropogenic gases and aerosols.

Regional Trends
To better examine the variability of the pollutants as well as regional trend differences, in Figure 17, we plotted the annual time series of all pollutants over six representative regions marked on Figure 1a. The thick black lines in Figure 17 show the four-year moving average smoothed time series. The percentage of change is also marked on Figure 17 with positive changes shown in red and negative in blue.
The changes of the time series in Figure 17 can mostly be inferred from the trend maps shown previously, but some detailed regional features are still observed. For example, although in general, the AOD has started to decrease since 2008, at NCP, the AOD increased until~2012 before a declining trend started. NO 2 showed similar behaviors for the NCP, YRD, and Sichuan Basin, where it first increased until~2012 and then experienced a sharp decline. Moreover, different to most of the country where overall negative AOD, NO 2 , and SO 2 trends in the last decade were found, over NW China, these parameters were steadily increasing, together with AAI. This may be partly associated with the China western development strategy, which increased the emission of anthropogenic gases and aerosols. As aerosol properties and pollution composition have different seasonal variability over different regions, we also examined the seasonal changes in detail. Figure 18 shows the seasonal time series for the six variables over the six representative regions. Although a glance of Figure 16 mostly indicates consistency with annual trends, more careful examination shows considerable differences in the seasonal trends. The difference in AOD trends tends to be larger than trace gases, which is reasonable, since AOD involves the complicated mixing of different species. Over the NCP, the start of the AOD decrease for the winter season (DJF) clearly lags behind the other three seasons, and the variability of the AOD time series appears more consistent with that of NO2 than those of SO2 and AAI. This suggests that changes in nitrate aerosols are likely to be responsible for the winter AOD trends in this region. The increase of AAI might also contribute to the lagged decline of winter aerosol loading and the near flat aerosol trend in the fall (SON). For the YRD, the strongest decline in AOD as large as 0.3/decade occurred in the summer season (JJA) after 2008; spring and fall also showed a moderate decrease, whereas the winter AOD time series was generally flat. Although SO2 and NO2 reduction can lead to a significant decrease of sulfate and nitrate aerosols, respectively, the increase of absorbing aerosols in the winter and fall as reflected in the positive AAI trends may compensate for this decrease. For the PRD region, summer and fall AOD exhibited a consistent decline from 2000 to 2018, whereas winter and spring AOD first increased until ~2008 and then decreased. This difference might be related to meteorology and biomass burning aerosols transported from Southeast Asia [60,70]. The Sichuan Basin showed consistent AOD changes except for winter, when a flat AOD time series was observed. Unlike the other three seasons, the winter season suffered from an increase in both AAI and AE, which indicated that the AOD decrease caused by SO2 and NO2 reduction might have been compensated by an increase in fine mode absorbing aerosols, mostly black and organic As aerosol properties and pollution composition have different seasonal variability over different regions, we also examined the seasonal changes in detail. Figure 18 shows the seasonal time series for the six variables over the six representative regions. Although a glance of Figure 16 mostly indicates consistency with annual trends, more careful examination shows considerable differences in the seasonal trends. The difference in AOD trends tends to be larger than trace gases, which is reasonable, since AOD involves the complicated mixing of different species. Over the NCP, the start of the AOD decrease for the winter season (DJF) clearly lags behind the other three seasons, and the variability of the AOD time series appears more consistent with that of NO 2 than those of SO 2 and AAI. This suggests that changes in nitrate aerosols are likely to be responsible for the winter AOD trends in this region. The increase of AAI might also contribute to the lagged decline of winter aerosol loading and the near flat aerosol trend in the fall (SON). For the YRD, the strongest decline in AOD as large as 0.3/decade occurred in the summer season (JJA) after 2008; spring and fall also showed a moderate decrease, whereas the winter AOD time series was generally flat. Although SO 2 and NO 2 reduction can lead to a significant decrease of sulfate and nitrate aerosols, respectively, the increase of absorbing aerosols in the winter and fall as reflected in the positive AAI trends may compensate for this decrease. For the PRD region, summer and fall AOD exhibited a consistent decline from 2000 to 2018, whereas winter and spring AOD first increased until~2008 and then decreased. This difference might be related to meteorology and biomass burning aerosols transported from Southeast Asia [60,70]. The Sichuan Basin showed consistent AOD changes except for winter, when a flat AOD time series was observed. Unlike the other three seasons, the winter season suffered from an increase in both AAI and AE, which indicated that the AOD decrease caused by SO 2 and NO 2 reduction might have been compensated by an increase in fine mode absorbing aerosols, mostly black and organic carbon. NW China is primarily dominated by dust aerosols, and the AOD values were clearly much higher in the dust-peak season of spring. However, the AOD for the non-peak seasons (i.e., summer, autumn, and winter) all increased at~0.03/decade while that for spring first increased rapidly at 0.1/decade until~2008, but then slightly decreased by~0.05 from 2008 to 2016. Meanwhile, the SO 2 and NO 2 concentrations also showed consistent increases, unlike East China. As discussed above, this might be related to the western development policy. An increase in anthropogenic aerosols should be responsible for the AOD increase during the non-peak seasons, whereas the dust change is most likely to be related to natural variability. The last region, NE China, is dominated by heavy industries and anthropogenic aerosols should be the major component. Unlike most East China regions, its AOD for all seasons except winter first showed clear decreases from 2000 to~2006 by~0.15, but remained flat afterward. Winter AOD showed weak and consistent upward trends (~0.04/decade). Additionally, unlike East China, no obvious reduction in SO 2 and NO 2 concentrations were found over this region, and NO 2 even increased for the winter and fall. Pollution control here might not have been implemented as strictly as East China, which led to the insignificant or even positive AOD trends. Finally, it should be noted that the changes in the loading of aerosols and trace gases can be associated with changes in both emissions and meteorology. In the future, we plan to carry out a series of detailed studies involving both data analysis and model simulations to attribute the trends.
Remote Sens. 2020, 12, x FOR PEER REVIEW 20 of 25 carbon. NW China is primarily dominated by dust aerosols, and the AOD values were clearly much higher in the dust-peak season of spring. However, the AOD for the non-peak seasons (i.e., summer, autumn, and winter) all increased at ~0.03/decade while that for spring first increased rapidly at 0.1/decade until ~2008, but then slightly decreased by ~0.05 from 2008 to 2016. Meanwhile, the SO2 and NO2 concentrations also showed consistent increases, unlike East China. As discussed above, this might be related to the western development policy. An increase in anthropogenic aerosols should be responsible for the AOD increase during the non-peak seasons, whereas the dust change is most likely to be related to natural variability.

Discussions
Our study is not the first to report the pollution trends in China. Previous studies have also noticed similar trends or change in AOD [56,71] or trace gases [20,72]. However, these studies mostly used a single set of satellite data or/and focused on a particular region. Our study was the first to conduct a comprehensive spatial analysis of pollution changes in China using multiple satellite

Discussions
Our study is not the first to report the pollution trends in China. Previous studies have also noticed similar trends or change in AOD [56,71] or trace gases [20,72]. However, these studies mostly used a single set of satellite data or/and focused on a particular region. Our study was the first to conduct a comprehensive spatial analysis of pollution changes in China using multiple satellite datasets. The consistency of different datasets not only increases the robustness of the results, but also provides valuable insights into the compositional changes. Compared to previous studies, our study offers a few new findings. (1) While previous studies have mostly focused on the most densely populated regions such as the NCP, YRD, and PRD, we showed that aerosol loading over NW and NE China has increased. This is most likely the result of increased dust, according to the positive AAI trends and negative AE trends. However, from the slight increase of SO 2 and NO 2 , we can infer that sulfate and nitrate aerosols may also play a role. (2) Many studies revealed the recent pollution reduction in East China, but the shift from positive to negative trends was different for different regions and for different pollutants. In particular, the pollution decline for the YRD and PRD started in~2008, whereas for the NCP and the Sichuan Basin, it started around~2012. A further comparison with SO 2 and NO 2 trends indicates that the YRD and PRD trends were more consistent with the SO 2 trends, while for the NCP and the Sichuan Basin, it was in phase with the NO 2 trends. (3) Previous study (e.g., [66]) about tropospheric O 3 trends have only used surface observations, therefore the analysis was focused on the east part because there were very few monitoring sites in West China. By using satellite data, we showed that there were also significant increases in tropospheric O 3 for West China, even over the Tibetan Plateau. The trends also exhibited distinct seasonal features. These findings highlight the usefulness and importance of satellites in environmental monitoring and research.
The cause of the observed trends, nonetheless, still needs investigation. The decrease of SO 2 , NO 2 , and AOD is likely to be the result of a series of emission reduction strategies taken by China. However, meteorology may also play an important role [73], especially concerning natural sources. In the future, we will be further quantifying different effects using chemical transport models. We will also examine the interannual variability and regional differences of pollution changes in detail.

Conclusions
In this study, we analyzed trends of major pollutants in China including aerosols, SO 2 , NO 2 , and tropospheric ozone using multiple satellite retrieval datasets. The major conclusions include: (1) Total aerosol loading in China, especially in the eastern and southern parts, has exhibited a significant decrease of~0.15-0.3/decade since 2008, after an almost decade-long increase. (2) This negative trend is accompanied by decreases of SO 2 and NO 2 during the same period, suggesting reductions in anthropogenic emissions. (3) NW and NE China, however, exhibit weak positive AOD trends at~0.1/decade, together with significant positive AAI and negative AE trends. These combined indicate that dust aerosol loading has likely increased and caused a total AOD increase over these two regions. (4) Multiple regression of AOD against AAI, SO 2 , and NO 2 show overall successful fitting and reasonable spatial distribution of the coefficients, further supporting the inferred compositional changes associated with the AOD trends. (5) Unlike aerosols, tropospheric ozone exhibits near uniform significant upward trends all over China, implying that ozone pollution control remains a challenging problem. ozone data were downloaded from https://acd-ext.gsfc.nasa.gov/Data_services/cloud_slice/#nd. The AERONET data for Beijing and XiangHe were provided by the AERONET website at https://aeronet.gsfc.nasa.gov/.

Conflicts of Interest:
The author declares no conflict of interests.