Joint Governance Regions and Major Prevention Periods of PM 2.5 Pollution in China Based on Wavelet Analysis and Concentration-Weighted Trajectory

: China has made some progress in controlling PM 2.5 (particulate matter with an aerodynamic diameter of ≤ 2.5 µ m) pollution, but there are still some key areas that need further strengthening. Considering that excessive prevention and control e ﬀ orts a ﬀ ect economic development, this paper combined an empirical orthogonal function, a continuous wavelet transform, and a concentration-weighted trajectory method to study joint regional governance during key pollution periods to provide suggestions for the e ﬃ cient control of PM 2.5 . The results from our panel of data of PM 2.5 in China from 2016 to 2018 could be decomposed into two modes. In the ﬁrst mode, the pollution center was in central Shaanxi Province, and the main eruption period was from November to January of the following year. As the center of this region, Xi’an should cooperate with the four cities in eastern Sichuan (Nanchong, Guangan, Bazhong, and Dazhou) to control PM 2.5 , since the eruption occurred in this area. Moreover, governance should last for at least two cycles, where one cycle is at least 23 days. The pollution center of the second mode was in the western part of Xinjiang. Therefore, after the prevention and control e ﬀ orts during the ﬁrst mode are completed, the regional city of Kashgar should continue to build a joint governance zone for PM 2.5 along the Tianshan mountains in the east, focusing on prevention and control over two cycles (where one cycle is 28 days).


Introduction
The occurrence of haze and severe PM 2.5 (particulate matter with an aerodynamic diameter of ≤2.5 µm) pollution has attracted broad attention around the world. Increased PM 2.5 concentrations lead to a deterioration in human health [1], visibility [2], the regional climate [3], and economic As for studies of polluted periods, previous researchers have mainly analyzed changes in PM 2.5 concentrations using statistical methods that include the year, quarter, or day to analyze variations in the concentration of PM 2.5 [7,9,15,21]. These researchers have used overall temporal information to analyze the evolutionary features of PM 2.5 based on a simple method, where temporal characteristics of PM 2.5 are evaluated with the assumption that the statistical properties of the time series do not vary over time. These types of studies can only uncover intense periods of high PM 2.5 , ignoring the pollution cycle. Because of that, pertinent suggestions for control cannot be given. Meanwhile, PM 2.5 concentrations have a complex cyclical variation with several short and long periods, which makes it difficult to analyze the temporal changes in PM 2.5 concentrations. The wavelet transform method is a feasible and effective method for studying the laws of variation of air pollution time-series indexes [22]. Temporal variations are expressed well, and mutated signals of the PM 2.5 level can be identified through wavelet analysis [23]. Using the wavelet method, Chen et al. [24] discovered multiscale features that are indicated in the temporal evolution of PM 2.5 . Huang et al. [25] used the wavelet transform method to obtain the characteristics of yearly changes and sudden changes in the PM 2.5 level. Liang et al. [26] utilized the wavelet approach to explore the potential association between PM 2.5 and influenza. In addition, the wavelet transform method is mainly used in PM 2.5 concentration predictions [27][28][29][30][31]. However, these studies have barely focused on the characteristics of periodic changes in the PM 2.5 concentration on a major scale.
The factors of time and space should both be considered in analyses of the period and location of heavy pollution. The empirical orthogonal function (EOF) method was originally introduced in meteorology as a method for extracting the dominant modes of spatial variability [32]: It has since been applied in climatological [33,34], hydrological [35,36], and geophysical studies [37]. This method can extract the main components of meteorological factors and temporal and spatial variations. This paper decomposes PM 2.5 concentrations and judges the pollution center through the PM 2.5 concentration distribution in different modes. Furthermore, it analyzes the counterpart time variation and cycle patterns, giving the foundation for an analysis of joint governance regions and control cycles.
Due to limitations with monitoring conditions, previous studies on the spatiotemporal distribution of PM 2.5 have mainly focused on the Yangtze River Delta [38], the Pearl River Delta [39,40], the Beijing-Tianjin-Hebei region [41], and other developed regions. Even if the scopes of these studies were national, pollutants in the western regions were characterized as missing or unpolluted due to missing data, which is not consistent with reality [7]. As the Chinese government has paid more attention to PM 2.5 issues, observation conditions across the country have improved. In 2013, there were only nearly 800 monitoring sites in the country [15], but as of December 2018, the number of sites exceeded 1600, covering the vast majority of the country. With these improvements in monitoring conditions, the pollution problem in the western region-which had not been studied before-has gradually come to light. Therefore, it is necessary to carry out a national PM 2.5 analysis of temporal and spatial variations according to existing conditions and to further expand previous research to create more accurate PM 2.5 studies in China.
First, this paper combines a spatiotemporal distribution of PM 2.5 with an analysis of the source of pollutants, studying the spatiotemporal distribution of PM 2.5 concentrations in China to identify seriously polluted areas and time coefficients. Moreover, we utilized a wavelet transform to analyze key pollution periods using time coefficients. Finally, a backward trajectory analysis and a concentration-weighted trajectory were utilized to study the joint prevention regions of seriously polluted areas at heavily polluted times, providing suggestions for the efficient control of PM 2.5 . The main contributions of this work are summarized as follows: • The study of the spatiotemporal distribution of PM 2.5 concentrations. We applied EOF decomposition to the spatiotemporal distribution of PM 2.5 concentrations, studying the time coefficients and vectors of PM 2.5 concentrations in different modes and analyzing the overall average state and local variation of PM 2.5 . Thus, regions and periods of heavy pollution could be determined. • An analysis of the duration of key prevention and control strategies. The time coefficient of PM 2.5 concentrations under different modalities was analyzed using a wavelet transform to judge the length of time of serious pollution periods so as to provide suggestions for the duration of prevention and control policies. • An investigation into the areas of joint protection and control. On the basis of the duration of key protection and control policies, PM 2.5 pollution in major cities in heavily polluted areas was analyzed using a backward trajectory analysis and a potential source analysis, with seriously polluted areas selected as the research object.
The rest of this paper is organized as follows ( Figure 1): Chapter 2 describes our data sources and research methods. Chapter 3.1 presents a descriptive statistical analysis of the annual variation in PM 2.5 from 2016 to 2018, using EOF to decompose panel data of PM 2.5 concentrations to identify the periods of and areas with heavy pollution in different modes. Chapter 3.2 uses a continuous wavelet transform to analyze the duration of key protection and control policies. Using the above-mentioned pollution areas and prevention periods, Chapter 3.3 analyzes the joint prevention and control areas of the research object from the point of view of air mass trajectory. Chapter 4 is the conclusion.
Sustainability 2020, 12, x FOR PEER REVIEW 4 of 23 • An investigation into the areas of joint protection and control. On the basis of the duration of key protection and control policies, PM2.5 pollution in major cities in heavily polluted areas was analyzed using a backward trajectory analysis and a potential source analysis, with seriously polluted areas selected as the research object. The rest of this paper is organized as follows ( Figure 1): Chapter 2 describes our data sources and research methods. Chapter 3.1 presents a descriptive statistical analysis of the annual variation in PM2.5 from 2016 to 2018, using EOF to decompose panel data of PM2.5 concentrations to identify the periods of and areas with heavy pollution in different modes. Chapter 3.2 uses a continuous wavelet transform to analyze the duration of key protection and control policies. Using the above-mentioned pollution areas and prevention periods, Chapter 3.3 analyzes the joint prevention and control areas of the research object from the point of view of air mass trajectory. Chapter 4 is the conclusion. The β-ray decay method and the tapered element oscillating microbalance (TEOM) method are commonly used as monitoring methods in China. Although these two methods are different in terms of measuring the concentration of PM2.5 and although the TEOM is known to have seasonally dependent biases [44], there is no evidence that China's ground monitoring data are not valid. Therefore, this paper used the (valid) ground monitoring data of PM2.5 concentrations released by official data sources, which are widely used in academic research [8,15,24], without considering errors in the measurements. This paper collected 24 h monitoring data from 361 cities in China, and daily and annual data were obtained through averaging. The positions of the selected cities are shown in Figure 2. An inverse distance weight algorithm [45] was employed to interpolate the PM2.5 concentration and a simulation of the spatial distribution of pollution. The results, combined with geographic data, are presented in ArcGIS 10.2.

Backward Trajectory Analysis Data
The meteorological data, which were provided by the NCEP (National Environmental Forecasting Center), were exploited in backward trajectory mode (GDAS1 (Global Data Assimilation System) data). The meteorological element field included temperature, air pressure, relative humidity, ground precipitation, horizontal and vertical wind speed, etc., from 2016 to 2018. GADS1 has the capability of calculating trajectory directly using vertical wind speed, which is an advantage over  . The β-ray decay method and the tapered element oscillating microbalance (TEOM) method are commonly used as monitoring methods in China. Although these two methods are different in terms of measuring the concentration of PM 2.5 and although the TEOM is known to have seasonally dependent biases [44], there is no evidence that China's ground monitoring data are not valid. Therefore, this paper used the (valid) ground monitoring data of PM 2.5 concentrations released by official data sources, which are widely used in academic research [8,15,24], without considering errors in the measurements. This paper collected 24 h monitoring data from 361 cities in China, and daily and annual data were obtained through averaging. The positions of the selected cities are shown in Figure 2. An inverse distance weight algorithm [45] was employed to interpolate the PM 2.5 concentration and a simulation of the spatial distribution of pollution. The results, combined with geographic data, are presented in ArcGIS 10.2.

Materials and Methods
Sustainability 2020, 12, x FOR PEER REVIEW 5 of 23 other methods, whose vertical wind speed is calculated indirectly by calculating the vertical integration of horizontal wind speed divergence [46].

Empirical Orthogonal Function
The empirical orthogonal function is a field analysis method widely used in the sphere of geosciences. Its principle is to decompose the spatiotemporal element field into several spatial basic modes and a linear combination of the time coefficient series, and then to objectively and quantitatively analyze the spatial structure and time changes of the element fields. The panel data of m observation points and n observations were expanded using the EOF and decomposed into the sum of the product of the orthogonal space matrix V and the orthogonal time matrix T: 11 1 Equation (1) is multiplied on the right to get Equation (2): where the superscript T represents the transpose of the matrix,  is a diagonal matrix composed of the eigenvalues of the matrix, and V is a matrix composed of the matrix eigenvectors. Therefore, the time coefficient can be obtained from Equation (3): (3) Figure 2. The distribution of the 361 cities used as data sources.

Backward Trajectory Analysis Data
The meteorological data, which were provided by the NCEP (National Environmental Forecasting Center), were exploited in backward trajectory mode (GDAS1 (Global Data Assimilation System) data). The meteorological element field included temperature, air pressure, relative humidity, ground precipitation, horizontal and vertical wind speed, etc., from 2016 to 2018. GADS1 has the capability of calculating trajectory directly using vertical wind speed, which is an advantage over other methods, whose vertical wind speed is calculated indirectly by calculating the vertical integration of horizontal wind speed divergence [46].

Empirical Orthogonal Function
The empirical orthogonal function is a field analysis method widely used in the sphere of geosciences. Its principle is to decompose the spatiotemporal element field into several spatial basic modes and a linear combination of the time coefficient series, and then to objectively and quantitatively analyze the spatial structure and time changes of the element fields. The panel data of m observation points and n observations were expanded using the EOF and decomposed into the sum of the product of the orthogonal space matrix V and the orthogonal time matrix T: Equation (1) is multiplied on the right to get Equation (2): where the superscript T represents the transpose of the matrix, Λ is a diagonal matrix composed of the eigenvalues of the matrix, and V is a matrix composed of the matrix eigenvectors. Therefore, the time coefficient can be obtained from Equation (3): North's Rule of Thumb assesses the uniqueness of EOF modes through assumptions of error of the eigenvalues. When adjacent eigenvalues match the condition λ j+1 − λ j ≥ λ j ( 2 n ) 1 2 , these two eigenvalues and their corresponding modes pass the sampling error test [47].
The variance contribution rate ρ of the eigenvalues and the cumulative variance contribution rate of the first p eigenvalues are calculated as follows: In the EOF, a set of eigenvalues, eigenvectors, and time coefficients represents a distribution mode. The first few eigenvectors passing the significance test represent the maximum distribution structure. The component of the eigenvector with the largest absolute value represents the intensity center. If the positive and negative signs in the eigenvector are consistent, the eigenvector reflects the features with the same change trend. If the component of a certain eigenvector is in a positive and negative phase distribution, then this eigenvector represents two opposite distribution types. The time coefficient represents the time variation characteristics of the spatial distribution form. When the time coefficient is positive, the year is consistent with the distribution form represented by the eigenvector, and vice versa. The larger the absolute value of the time coefficient is, the more significant the distribution form. The variance contribution rate ρ reflects the degree to which a mode explains the whole. The higher its cumulative value is, the more accurately the selected mode describes the overall situation.

Continuous Wavelet Transform
A continuous wavelet transform can clearly reveal the period of change hidden in nonstationary time series and can reflect its change trend in different time scales. A continuous wavelet transform can also extract multiple wave periods from the wave sequence at different scales to reflect its changing trend, which is suitable for signal feature extraction [48]. The principle is to shift the mother wavelet function ψ(t) after the translation b, and then to make an inner product with the signal f (t) for it to be analyzed at different scales a, as follows: where b = 1, 2, . . . , N. N is the number of datapoints (1096 in this paper). According to the literature [49], for nonorthogonal wavelet analyses, one can use an arbitrary set of scales a to build a more complete picture. It is convenient to write the scales as fractional powers of two: where s 0 is the smallest resolvable scale and J determines the largest scale. Here, δt represents the sampling intervals of data in this paper, and s 0 should be chosen so that the equivalent Fourier period is approximately 2δt. For the Morlet wavelet, a δ j of about 0.5 is the largest value that can still give an adequate sampling scale. Parameter values from the literature [49] were used here: δ j = 0.125. Therefore, the smallest and largest resolvable scales could be calculated as 2 and 512, respectively. Note that longer scales correspond to the most stretched-out wavelets. The more stretched-out the wavelet is, the longer the portion of the signal to which it is being compared (and the coarser the signal features measured by the wavelet coefficients) will be. In addition, a large scale also means a long period. The time between the adjacent maximum and minimum of wavelet coefficients is almost 7 months when the scale a is over 300 d (with similar data conditions) [24]. Because this article focuses on short-term trends in PM 2.5 concentration changes for efficient governance, the largest scale was set at 64, the same as in the literature [49]. The complex Morlet wavelet, which is a single-frequency complex sinusoidal function, which has symmetry, nonorthogonality, and an imaginary part [24], can simultaneously preserve the amplitude and phase information of the sequence signal of PM 2.5 concentrations. Therefore, a complex Morlet wavelet was selected as the mother wavelet. The complex Morlet wavelet function is as follows: where t represents time and f b represents the bandwidth controlling attenuation in the time domain and the corresponding bandwidth in the frequency domain. Here, f b is the reciprocal of variance in the frequency domain. An increase in f b will lead to the wavelet energy being concentrated around the center frequency and will slow down the attenuation speed in the time domain. On the other hand, a decrease in f b will accelerate the decay rate in the time domain and reduce the energy of the frequency domain. Here, f c denotes the center frequency and affects the frequency value when the time domain is converted into a frequency domain [24,50]. These two parameters can be adjusted to obtain appropriate time-frequency resolutions. In this paper, we set the f b as 1 and the f c as 1.5, in accordance with parameter settings from the literature [24] and parameter optimization from the literature [51].
Using the MatlabR2018a software platform, the wavemenu toolbox was selected to calculate the wavelet coefficients. Meanwhile, the wavelet variance of the EOF time coefficient was analyzed to find the main scale. Then, the coefficient modulus of the wavelet coefficients and the contour plots of the real part were plotted to find the main oscillation periods and periods of time series at major scales. The periodic variation of the real part of the wavelet coefficient on the main scale was the most significant period of the time series.
A continuous wavelet transform can expand one-dimensional signals in both the time and frequency directions to analyze the time-frequency structure of the data in detail and extract valuable information. A continuous wavelet transform can not only provide the relative contribution of different scales of the time series, but can also indicate changes in different scales, so it is very helpful for the study of time series.

Backward Trajectory Analysis and Concentration-Weighted Trajectory
Here, we convey the calculation of the air trajectory first. The backward trajectory mode adopts the Hybrid Single-Particle Lagrangian Integrated Trajectory Model (HYSPLIT) developed by the NOAA (National Oceanic and Atmospheric Administration) [52,53]. The HYSPLIT model uses gridding meteorological data to respond to emergencies in the atmosphere to diagnose issues and analyze the climate. This mode is Lagrangian, as is Euler's mixed diffusion mode: The processes of advection and diffusion are calculated using the Lagrangian method, and the concentration is calculated using the Euler Method [54]. HYSPLIT is considerably detailed in terms of the transportation, diffusion, and settling of pollutants, and its highest simulation accuracy can last hours. Thus, HYSPLIT is widely used in analyzing the source of pollutants and determining transmission and diffusion [55].
MeteoinfoMap is a geographic information system application that analyzes and visualizes multiple meteorological data formats [56]. Its plug-in, TrajStat, can use the backward trajectory analysis software adopted by the HYSPLIT Lagrangian diffusion module [57]. On the basis of the results from the EOF analysis, we selected cities with high pollution levels as the starting point and used the length of time of severe pollution as the pushback time. The air mass movement path at 500 m of altitude (the wind field at 500 m of altitude) can reflect the characteristics of the average flow field in the boundary layer [16,17,19]. Therefore, the height of the simulation was selected as 500 m, and we calculated the 48 h backward trajectory every 2 h.
After that, a grid was established based on the area covered by the atmospheric trajectory, with a resolution of 0.25 • × 0.25 • . The concentration-weighted trajectory (CWT) was used to analyze the source of pollution. The CWT is a mixed-trajectory receptor model that combines meteorological trajectory nodes (residence time) and pollutant concentrations to trace their contributions to the pollution of a recipient site [58]. After the study area was gridded, the CWT value of Grid (i, j) was as follows: where τ ijl is the number of pollution trajectory nodes in the grid (i, j) in the area, and C l is the pollutant concentration of the trajectory. The higher the grid CWT value is, the greater the probability that the pollution trajectory comes from that grid point. When there are fewer trajectories within the grid, the residence time is shorter, so the x value is higher and there is greater uncertainty. When n ij is less than three times its average value, we use the following weight function to reduce the uncertainty [58][59][60]: 1.00 3n ave < n ij 0.70 n ave < n ij ≤ 3n ave 0.42 0.4n ave < n ij ≤ n ave 0.05 n ij ≤ 0.4n ave

Descriptive Statistical Analyses
According to the PM 2.5 concentration data from 2016 to 2018, the annual mean of the PM 2.5 concentration across the country had a downward trend. The national average PM 2.5 concentration value was reduced by a total of 6.753 µg·m -3 , representing 14.6% of the former value. The mean value decreased from 46.15 µg·m -3 in 2016 to 43.76 µg·m -3 in 2017, and then declined to 39.39 µg·m -3 . As is shown in Figure 3, except for the North China Plain, southern Sichuan, and part of western Xinjiang, the three-year average PM 2.5 concentration value in most areas was less than 55 µ·m _3 , which indicates good-quality air. In 2016, heavily polluted areas (in terms of PM 2.5 pollution) were mainly concentrated in western Xinjiang, the Beijing-Tianjin-Hebei region, Shanxi, Henan, central Hubei, central Shanxi, and west Shandong. In 2017, the concentration of PM 2.5 in the Shanxi, Anhui, Guangxi, Guangdong, Hainan, and Northern Xinjiang areas increased. At the same time, the PM 2.5 concentration in western Xinjiang dropped significantly. However, due to the high PM 2.5 concentration in western Xinjiang, although the decline was large, the air quality in Xinjiang was still worrying. Although PM 2.5 concentrations in most parts of the country were under control in 2018, PM 2.5 concentrations in western Xinjiang and Tibet continued to rebound slightly, and PM 2.5 concentrations also increased in western Yunnan and northern Gansu, indicating that in these areas, regional PM 2.5 prevention and control policies were not effective. These results differ from previous studies of the PM 2.5 distribution in China (by Hu Line [41]).
Sustainability 2020, 12, x FOR PEER REVIEW 9 of 23 problem in some regions. PM2.5 concentrations did not always continue to decrease over the three years. Therefore, PM2.5 prevention and control work should focus on these problem areas by strengthening control policies further, improving governance policies, and making PM2.5 governance more effective.

Heavily Polluted Areas and Periods
In order to further understand the spatial and temporal distribution of PM2.5 pollution and to accurately determine the key areas and periods of PM2.5 governance, this paper used EOF analysis. We employed a standardized matrix of monthly panel data of PM2.5 concentrations in 361 cities from 2016 to 2018. The characteristic vector of EOF expansion demonstrated the spatial distribution of the PM2.5 concentration in each province. The maximum value center of all eigenvector fields indicated the region that was the most sensitive to PM2.5 concentration changes. The eigenvectors of PM2.5 concentrations in China's provinces converged quickly, with eight eigenvalues passing North's Rule of Thumb. The variance contribution rates of each mode were 41.70%, 10.77%, 6.56%, 4.72%, 3.63%, 2.93%, 2.12%, and 1.78%. The contribution rates of the eigenvalue variance illustrated that the EOF analysis results were ideal. Compared to the first two modes, the third mode and the subsequent These research results should alarm the Chinese government: Although PM 2.5 has been under control in most regions of the country under current policies, PM 2.5 pollution remains a major problem in some regions. PM 2.5 concentrations did not always continue to decrease over the three years. Therefore, PM 2.5 prevention and control work should focus on these problem areas by strengthening control policies further, improving governance policies, and making PM 2.5 governance more effective.

Heavily Polluted Areas and Periods
In order to further understand the spatial and temporal distribution of PM 2.5 pollution and to accurately determine the key areas and periods of PM 2.5 governance, this paper used EOF analysis. We employed a standardized matrix of monthly panel data of PM 2.5 concentrations in 361 cities from 2016 to 2018. The characteristic vector of EOF expansion demonstrated the spatial distribution of the PM 2.5 concentration in each province. The maximum value center of all eigenvector fields indicated the region that was the most sensitive to PM 2.5 concentration changes. The eigenvectors of PM 2.5 concentrations in China's provinces converged quickly, with eight eigenvalues passing North's Rule of Thumb. The variance contribution rates of each mode were 41.70%, 10.77%, 6.56%, 4.72%, 3.63%, 2.93%, 2.12%, and 1.78%. The contribution rates of the eigenvalue variance illustrated that the EOF analysis results were ideal. Compared to the first two modes, the third mode and the subsequent modes took up a small proportion of the EOF decomposition. Therefore, the first two eigenvectors could be used to represent the spatiotemporal structure of PM 2.5 concentrations.
The variance contribution of the first mode of EOF indicated the average state of the PM 2.5 concentration from 2016 to 2018, representing the distribution field of PM 2.5 concentrations in China. In Figure 4, it can be seen that the first eigenvector had a consistently positive value, indicating that the spatial variation trend of PM 2.5 concentrations in the first mode was synchronous. The distribution pattern of the eigenvector in the first mode was close to the mean distribution (Figure 3a). It can be concluded that the first mode was the average state of PM 2.5 concentration. The highest value area was located in the middle of Shaanxi Province. The eigenvalues of North China, central Sichuan, and parts of Xinjiang were also large. This indicates that the PM 2.5 concentration fluctuated greatly in these areas. When the time coefficient was positive, PM 2.5 pollution was more severe. In other areas, the eigenvalue was close to zero regardless of changes in the time coefficient, so the PM 2.5 concentration remained at a low level.
Sustainability 2020, 12, x FOR PEER REVIEW 10 of 23 modes took up a small proportion of the EOF decomposition. Therefore, the first two eigenvectors could be used to represent the spatiotemporal structure of PM2.5 concentrations.
The variance contribution of the first mode of EOF indicated the average state of the PM2.5 concentration from 2016 to 2018, representing the distribution field of PM2.5 concentrations in China. In Figure 4, it can be seen that the first eigenvector had a consistently positive value, indicating that the spatial variation trend of PM2.5 concentrations in the first mode was synchronous. The distribution pattern of the eigenvector in the first mode was close to the mean distribution (Figure 3a). It can be concluded that the first mode was the average state of PM2.5 concentration. The highest value area was located in the middle of Shaanxi Province. The eigenvalues of North China, central Sichuan, and parts of Xinjiang were also large. This indicates that the PM2.5 concentration fluctuated greatly in these areas. When the time coefficient was positive, PM2.5 pollution was more severe. In other areas, the eigenvalue was close to zero regardless of changes in the time coefficient, so the PM2.5 concentration remained at a low level. From the distribution of the time coefficients in the first mode ( Figure 5), it can be concluded that the PM2.5 concentration change in the first mode had obvious seasonal characteristics. In the first and fourth quarters of each year, most PM2.5 time coefficients were positive, and the maximum value each year appeared around December, which means that the high-value area of the eigenvector in the first and fourth quarters represented the PM2.5 pollution problem getting worse and the pollution being the most serious from December to January. In the second and third quarters, the time coefficients were negative, so the PM2.5 concentration gradually decreased. The dashed line in the From the distribution of the time coefficients in the first mode ( Figure 5), it can be concluded that the PM 2.5 concentration change in the first mode had obvious seasonal characteristics. In the first and fourth quarters of each year, most PM 2.5 time coefficients were positive, and the maximum value each year appeared around December, which means that the high-value area of the eigenvector in the first and fourth quarters represented the PM 2.5 pollution problem getting worse and the pollution being the most serious from December to January. In the second and third quarters, the time coefficients were negative, so the PM 2.5 concentration gradually decreased. The dashed line in the figure is the linear regression of the time coefficients of the first mode, which had a downward trend, reflecting that in the long term, the pollution problems in the areas where the PM 2.5 concentrations changed after the first mode were alleviated.
Sustainability 2020, 12, x FOR PEER REVIEW 11 of 23 figure is the linear regression of the time coefficients of the first mode, which had a downward trend, reflecting that in the long term, the pollution problems in the areas where the PM2.5 concentrations changed after the first mode were alleviated. According to the above analysis, in the first mode, the main pollution areas were Shaanxi, North China, Sichuan, and Xinjiang, and central Shaanxi. The main pollution periods were the first and fourth quarters of each year, and pollution was the most serious in December. Governments should focus on the above areas and pollution periods. The second EOF mode reflected the spatial variation of regional PM2.5 concentration differences; the spatiotemporal distributions are shown in Figure 6. Most of the eigenvalues exhibited in Figure  6 were close to zero, which meant that these regions contributed little to local pollution characteristics. However, the eigenvalues in Western Xinjiang had large negative values, suggesting dramatic changes in the PM2.5 concentration.  According to the above analysis, in the first mode, the main pollution areas were Shaanxi, North China, Sichuan, and Xinjiang, and central Shaanxi. The main pollution periods were the first and fourth quarters of each year, and pollution was the most serious in December. Governments should focus on the above areas and pollution periods.
The second EOF mode reflected the spatial variation of regional PM 2.5 concentration differences; the spatiotemporal distributions are shown in Figure 6. Most of the eigenvalues exhibited in Figure 6 were close to zero, which meant that these regions contributed little to local pollution characteristics. However, the eigenvalues in Western Xinjiang had large negative values, suggesting dramatic changes in the PM 2.5 concentration.
It can be seen from Figure 7 that the time coefficients had large negative values from February to April each year, reaching a maximum on April 3, 2018. When looking at the eigenvectors of western Xinjiang in the second mode, in areas where the eigenvector is negative, negative time coefficients reflect a positive increase in the PM 2.5 concentration. The greater the absolute values of the time coefficients and eigenvalues are, the more severe the pollution is. It can be seen that the PM 2.5 pollution in the second mode was concentrated from February to April each year, and the pollution problem was extremely serious, especially in April (when there was an eruption of severe pollution). PM 2.5 governance in western Xinjiang should be strengthened during this period; the pollution problem in April is especially worthy of more attention.
The maximum value is marked as a square.
The second EOF mode reflected the spatial variation of regional PM2.5 concentration differences; the spatiotemporal distributions are shown in Figure 6. Most of the eigenvalues exhibited in Figure  6 were close to zero, which meant that these regions contributed little to local pollution characteristics. However, the eigenvalues in Western Xinjiang had large negative values, suggesting dramatic changes in the PM2.5 concentration.  It can be seen from Figure 7 that the time coefficients had large negative values from February to April each year, reaching a maximum on April 3, 2018. When looking at the eigenvectors of western Xinjiang in the second mode, in areas where the eigenvector is negative, negative time coefficients reflect a positive increase in the PM2.5 concentration. The greater the absolute values of the time coefficients and eigenvalues are, the more severe the pollution is. It can be seen that the PM2.5 pollution in the second mode was concentrated from February to April each year, and the pollution problem was extremely serious, especially in April (when there was an eruption of severe pollution). PM2.5 governance in western Xinjiang should be strengthened during this period; the pollution problem in April is especially worthy of more attention.

Major Prevention Period
The PM2.5 pollution periods in different modes were analyzed above. Local governments should focus primarily on severely polluted periods, especially in the western region of China, to strictly control PM2.5 and further strengthen prevention and control efforts. However, considering the economic costs brought on by strengthening prevention and control efforts, the periods that are the most polluted should be of higher priority to be the most efficient. Since changes in PM2.5 levels also have high-intensity periods, we suggest carrying out at least one complete control cycle for polluted periods. Prevention during heavily polluted times is key to helping improve control efficiency. Therefore, a wavelet transform was performed on the time coefficients in the two modes to further study what the best prevention time periods were.
The variance of wavelet transform coefficients can determine the main scale of the wavelet transform, which is plotted in Figure 8. The line chart represents the relative intensity of each scale in a time series. The scale at the corresponding peak is called the main time scale of the wavelet transform, and the periodic oscillation at this scale is the most significant. It is clear from Figure 8

Major Prevention Period
The PM 2.5 pollution periods in different modes were analyzed above. Local governments should focus primarily on severely polluted periods, especially in the western region of China, to strictly control PM 2.5 and further strengthen prevention and control efforts. However, considering the economic costs brought on by strengthening prevention and control efforts, the periods that are the most polluted should be of higher priority to be the most efficient. Since changes in PM 2.5 levels also have high-intensity periods, we suggest carrying out at least one complete control cycle for polluted periods. Prevention during heavily polluted times is key to helping improve control efficiency. Therefore, a wavelet transform was performed on the time coefficients in the two modes to further study what the best prevention time periods were.
The variance of wavelet transform coefficients can determine the main scale of the wavelet transform, which is plotted in Figure 8. The line chart represents the relative intensity of each scale in a time series. The scale at the corresponding peak is called the main time scale of the wavelet transform, and the periodic oscillation at this scale is the most significant. It is clear from Figure 8 that the wavelet variance of the time coefficients in the first mode had two distinct peaks corresponding to the time scales of 37 and 23 d. The periodic oscillation of the time scale around 37 d was the largest, which could mean that this is the main scale. In the second mode, the time scale of about 44 d had the largest variance. On the basis of the wavelet coefficient modulus and variance, the main scale of the time coefficient of changes in PM2.5 concentration in the two modes (and their corresponding time range) was found. The main scale was used in the frequency domain analysis to determine the time domain change at different frequencies, which has no practical meaning. Because we wanted to focus on the most polluted periods from 2016 to 2018, we selected the scale with the largest variance of wavelet coefficients and the time range of the corresponding maximum value of the wavelet coefficient modulus. In practical applications, other typical scales or periods can also be selected for analysis according to actual needs, or all periods corresponding to the main scales can be selected for a comprehensive analysis. This method is very flexible. A contour plot of the real part of the wavelet transform coefficient can reflect the change in the data at different transform scales. The ordinate of the densely arranged region at the center of the positive and negative values is the characteristic time scale of the wave sequence, and the abscissa is the periodic time frame of changes. In Figure 10, there are four red and blue local oscillation regions in the first mode and two local oscillation regions in the second mode. The two types of oscillations at different scales reflect different periodic changes. Each time the regions alternate between red and blue is a pollution period. The above analysis was carried out by selecting areas located in the main time scale and within the strong oscillating periods. The first mode (Figure 10a, the 37 d scale) experienced two pronounced concentration bursts between December 2017 and January 2018. In the second mode (Figure 10b), the cold-and warm-colored alternation regions were mainly concentrated in the period between February and April in 2016, with two cycles of strength and weakness on a 44 On the basis of the wavelet coefficient modulus and variance, the main scale of the time coefficient of changes in PM2.5 concentration in the two modes (and their corresponding time range) was found. The main scale was used in the frequency domain analysis to determine the time domain change at different frequencies, which has no practical meaning. Because we wanted to focus on the most polluted periods from 2016 to 2018, we selected the scale with the largest variance of wavelet coefficients and the time range of the corresponding maximum value of the wavelet coefficient modulus. In practical applications, other typical scales or periods can also be selected for analysis according to actual needs, or all periods corresponding to the main scales can be selected for a comprehensive analysis. This method is very flexible. A contour plot of the real part of the wavelet transform coefficient can reflect the change in the data at different transform scales. The ordinate of the densely arranged region at the center of the positive and negative values is the characteristic time scale of the wave sequence, and the abscissa is the periodic time frame of changes. In Figure 10, there are four red and blue local oscillation regions in the first mode and two local oscillation regions in the second mode. The two types of oscillations at different scales reflect different periodic changes. Each time the regions alternate between red and blue is a pollution period. The above analysis was carried out by selecting areas located in the main time scale and within the strong oscillating periods. The first mode (Figure 10a, the 37 d scale) experienced two pronounced concentration bursts between December 2017 and January 2018. In the second mode (Figure 10b), the cold-and warm-colored alternation regions were mainly concentrated in the period between February and April in 2016, with two cycles of strength and weakness on a 44 On the basis of the wavelet coefficient modulus and variance, the main scale of the time coefficient of changes in PM 2.5 concentration in the two modes (and their corresponding time range) was found. The main scale was used in the frequency domain analysis to determine the time domain change at different frequencies, which has no practical meaning. Because we wanted to focus on the most polluted periods from 2016 to 2018, we selected the scale with the largest variance of wavelet coefficients and the time range of the corresponding maximum value of the wavelet coefficient modulus. In practical applications, other typical scales or periods can also be selected for analysis according to actual needs, or all periods corresponding to the main scales can be selected for a comprehensive analysis. This method is very flexible.
A contour plot of the real part of the wavelet transform coefficient can reflect the change in the data at different transform scales. The ordinate of the densely arranged region at the center of the positive and negative values is the characteristic time scale of the wave sequence, and the abscissa is the periodic time frame of changes. In Figure 10, there are four red and blue local oscillation regions in the first mode and two local oscillation regions in the second mode. The two types of oscillations at different scales reflect different periodic changes. Each time the regions alternate between red and blue is a pollution period. The above analysis was carried out by selecting areas located in the main time scale and within the strong oscillating periods. The first mode (Figure 10a, the 37 d scale) experienced two pronounced concentration bursts between December 2017 and January 2018. In the second mode (Figure 10b), the cold-and warm-colored alternation regions were mainly concentrated in the period between February and April in 2016, with two cycles of strength and weakness on a 44 d scale. It is worth noting that there was a slight increase in the concentration of PM 2.5 before the eruption of the two modes, which indicates the beginning of a serious eruption of PM 2.5 pollution, an important early warning sign for PM 2.5 prevention. After the second eruption, the PM 2.5 concentration decreased to normal levels, and we defined the end of the second eruption as the end time of major pollution.  The above analysis determined the main scales, strong oscillation periods, and cyclic changes in the corresponding time coefficients. For the purpose of determining the PM2.5 pollution cycle in different modes, the real part of the wavelet transform coefficients of different modes (with their main scales) was drawn. A typical polluted period was calculated using the intervals between the PM2.5 concentration peaks. In Figure 11, the peak is the eruption time, the start is the short eruption time before the intense eruption of pollution, and the end is the ending time of the major pollution period. With these calculations, we could find out that the typical pollution time under the first mode was between 27 November 2017 and 26 January 2018, the outbreak interval was 23 days, and there were two intense pollution eruptions of PM2.5. Pollution in the second mode occurred between 27 January 2016 and 31 March 2016, the outbreak interval was 28 days, and there were two intense pollution eruptions of PM2.5 here as well. The above analysis determined the main scales, strong oscillation periods, and cyclic changes in the corresponding time coefficients. For the purpose of determining the PM 2.5 pollution cycle in different modes, the real part of the wavelet transform coefficients of different modes (with their main scales) was drawn. A typical polluted period was calculated using the intervals between the PM 2.5 concentration peaks. In Figure 11, the peak is the eruption time, the start is the short eruption time before the intense eruption of pollution, and the end is the ending time of the major pollution period. With these calculations, we could find out that the typical pollution time under the first mode was between 27 November 2017 and 26 January 2018, the outbreak interval was 23 days, and there were two intense pollution eruptions of PM 2.5 . Pollution in the second mode occurred between 27 January 2016 and 31 March 2016, the outbreak interval was 28 days, and there were two intense pollution eruptions of PM 2.5 here as well.
Autocorrelation was used to verify the wavelet transform results. The autocorrelation results from the first and second modes' time coefficients are shown in Figures A1 and A2, respectively. According to the literature [61,62], the time differences between peaks can be calculated-we did that here to find that there existed a long period of 22 days, a result very close to the wavelet transform results from the first mode. Though the autocorrelation of the second mode was weaker than that of the first mode, a long period of 25 days was found for the time coefficient of the second mode, which was only a three-day difference compared to the wavelet analysis results. The autocorrelation results certified that the continuous wavelet transform method could find periods of time series from weak autocorrelation and can be used reliably in PM 2.5 concentration research. different modes, the real part of the wavelet transform coefficients of different modes (with their main scales) was drawn. A typical polluted period was calculated using the intervals between the PM2.5 concentration peaks. In Figure 11, the peak is the eruption time, the start is the short eruption time before the intense eruption of pollution, and the end is the ending time of the major pollution period. With these calculations, we could find out that the typical pollution time under the first mode was between 27 November 2017 and 26 January 2018, the outbreak interval was 23 days, and there were two intense pollution eruptions of PM2.5. Pollution in the second mode occurred between 27 January 2016 and 31 March 2016, the outbreak interval was 28 days, and there were two intense pollution eruptions of PM2.5 here as well. Autocorrelation was used to verify the wavelet transform results. The autocorrelation results from the first and second modes' time coefficients are shown in Figures A1 and A2, respectively. According to the literature [61,62], the time differences between peaks can be calculated-we did that here to find that there existed a long period of 22 days, a result very close to the wavelet transform results from the first mode. Though the autocorrelation of the second mode was weaker than that of the first mode, a long period of 25 days was found for the time coefficient of the second mode, which On the basis of the above results, Shaanxi, northern China, Sichuan, and some parts of Xinjiang (whose PM 2.5 performances were poor during the first mode) should commence their governance between November and January of the following year (when PM 2.5 pollution is severe). The governance in these regions should last for two cycles after slight PM 2.5 pollution occurs for the first time, and one cycle should last for no less than 23 days. Governance in western Xinjiang should be strengthened in terms of prevention and control (in the first mode) until March. There should be at least two cycles during this period, and one cycle should last for 28 days.

Joint Governance Region
On the basis of the analysis above, we selected the cities of Xi'an and Kashgar for analysis. Both cities were in the most polluted areas in both modes. Taking the most severely polluted period as the calculation time, we used a backward trajectory of 48 h to analyze the major joint-governance cities. The backward trajectory analysis and concentration-weighted trajectory conditions are shown in Table 1. Here, the atmospheric trajectory image and CWT values reflect the effect of grid points on the PM 2.5 levels of the analyzed city. The larger the grid value is, the more serious the impact of PM 2.5 pollution on the city. In order to better identify major pollution sources, only high-value CWT grids in China are plotted, which are shown in Figures 12 and 13 Here, the atmospheric trajectory image and CWT values reflect the effect of grid points on the PM2.5 levels of the analyzed city. The larger the grid value is, the more serious the impact of PM2.5 pollution on the city. In order to better identify major pollution sources, only high-value CWT grids in China are plotted, which are shown in Figure 12 and Figure 13. The CWT diagram of all grids is shown in Figures A3 and A4.  In Figure 12, it can be seen that Nanchong, Guangan, Bazhong, and Dazhou in the eastern Sichuan Province had the greatest impact on PM2.5 pollution in Xi'an. During the most polluted winter, PM2.5 in eastern Sichuan passed through Ankang, Hanzhong, and Shangluo in southern Shaanxi Province and reached Xi'an, exacerbating PM2.5 pollution. In addition, Weinan, Yanan, and Yuncheng in the northeast were also sources of PM2.5 pollution in Xi'an. Therefore, prevention and control strategies in Xi'an should not be limited to local governance, and joint governance should be carried out between the above cities to control pollution more effectively.
The sources of severe pollution in Kashgar are shown in Figure 13. The PM2.5 concentration in Kashgar was most influenced by the eastern region, which had heavy pollution. The grids of these trajectories were distributed along the Tianshan Mountains, starting from Bayingol and reaching Kashgar via Aksu, Aral, Tumshuk, and North Kashgar. Therefore, Kashgar should create prevention and control strategies for the urban agglomerations along the Tianshan Mountains.

Conclusions
This paper analyzed joint governance regions and key governing cycles in China through the EOF and backward trajectory analysis. The EOF decomposed PM2.5 concentration panel data into space vectors and time coefficients. The time coefficients reflected the changing trends in pollution distribution. After transforming the time coefficients using wavelets, we could obtain the pollution periods of the main time scales and could determine key governance periods using the real parts of the wavelet coefficients. A space vector indicated the distribution of PM2.5 pollution, which helped identify areas with serious pollution problems. Then, we selected the main cities in these areas and analyzed their sources of pollution during key governance periods to determine possible joint prevention areas. The conclusions can be summarized as follows: • From a nationwide point of view, the overall PM2.5 pollution in China improved from 2016 to 2018, but the improvement was limited. PM2.5 remains a serious problem in Xinjiang and North China. This finding is different from previous studies that have claimed that China's PM2.5 In Figure 12, it can be seen that Nanchong, Guangan, Bazhong, and Dazhou in the eastern Sichuan Province had the greatest impact on PM 2.5 pollution in Xi'an. During the most polluted winter, PM 2.5 in eastern Sichuan passed through Ankang, Hanzhong, and Shangluo in southern Shaanxi Province and reached Xi'an, exacerbating PM 2.5 pollution. In addition, Weinan, Yanan, and Yuncheng in the northeast were also sources of PM 2.5 pollution in Xi'an. Therefore, prevention and control strategies in Xi'an should not be limited to local governance, and joint governance should be carried out between the above cities to control pollution more effectively.
The sources of severe pollution in Kashgar are shown in Figure 13. The PM 2.5 concentration in Kashgar was most influenced by the eastern region, which had heavy pollution. The grids of these trajectories were distributed along the Tianshan Mountains, starting from Bayingol and reaching Kashgar via Aksu, Aral, Tumshuk, and North Kashgar. Therefore, Kashgar should create prevention and control strategies for the urban agglomerations along the Tianshan Mountains.

Conclusions
This paper analyzed joint governance regions and key governing cycles in China through the EOF and backward trajectory analysis. The EOF decomposed PM 2.5 concentration panel data into space vectors and time coefficients. The time coefficients reflected the changing trends in pollution distribution. After transforming the time coefficients using wavelets, we could obtain the pollution periods of the main time scales and could determine key governance periods using the real parts of the wavelet coefficients. A space vector indicated the distribution of PM 2.5 pollution, which helped identify areas with serious pollution problems. Then, we selected the main cities in these areas and analyzed their sources of pollution during key governance periods to determine possible joint prevention areas. The conclusions can be summarized as follows: • From a nationwide point of view, the overall PM 2.5 pollution in China improved from 2016 to 2018, but the improvement was limited. PM 2.5 remains a serious problem in Xinjiang and North China. This finding is different from previous studies that have claimed that China's PM 2.5 pollution is strong in the east and weak in the west. National PM 2.5 concentration panel data were decomposed using the EOF, resulting in two modes. The first mode reflected the average state of PM 2.5 pollution in China. Seriously polluted areas included northern China, western Sichuan, and parts of Xinjiang. The most polluted areas were in central Shaanxi Province. Pollution in the first mode had significant seasonal characteristics, indicating that pollution in the first and fourth quarters was serious and that the pollution degree decreased in the second and third quarters. The second mode reflected the local pollution characteristics of the Xinjiang region, indicating that there was severe pollution from February to April.

•
Major prevention periods during the two EOF modes were studied using a continuous wavelet transform. This showed that the first mode had a typical oscillation, with a scale of 37 d, from November to January of the following year, and the main oscillation scale of the second mode was 44 d, from February 2016 to April 2016. In the areas where PM 2.5 pollution was in the first mode, prevention and control strategies should be carried out from November to January of the following year. After a small eruption in pollution occurs, prevention and control should be implemented for a period of no less than 23 days. After the end of the control cycle in the first mode, management and control of PM 2.5 pollution in the second mode should be strengthened for at least two cycles until March, where one cycle is 28 days. • This paper took Xi'an and Kashgar as examples to analyze joint governance regions based on pollution trajectories. In addition to local control, PM 2.5 control strategies during winter in Xi'an should include joint control with Nanchong, Guangan, Bazhong, and Dazhou in eastern Sichuan. PM 2.5 control in winter in Xi'an should also be coordinated with Ankang, Hanzhong, and Shangluo in southern Shaanxi to alleviate the problem of pollution sources in the south. Meanwhile, Xi'an should work with Weinan, Yanan, and Yuncheng to solve the problem of pollution sources in the northwest. For Kashgar, it is necessary to establish a PM 2.5 joint governance area along the Tianshan Mountains to the east, with a focus on joint governance with the northern cities of Kashgararea, Tumshuk, Aksuarea, and Aral in the east.
The method proposed in this paper can be used to analyze the joint governance regions and periods of PM 2.5 pollution in key cities across the entire country, and it can also be used in smaller areas, such as provinces or cities. The method is simple, feasible, and has universality. However, due to the many factors affecting the concentration of PM 2.5 and the limited space of this paper, the important influencing factors of PM 2.5 pollution prevention and control were not analyzed, which should be the focus of further research.        Figure A4. CWTs of all grids in the second mode. The area of analysis of this article is within the purple rectangle.