Analysis of Synergistic Changes in PM_2.5 and O₃ Concentrations Based on Structural Equation Model Study

Abstract

Given the increasing importance of effectively identifying synergistic changes between PM_2.5 and O₃ and comprehensively analyzing their impact on air quality management in China, we employ the Sen+Mann–Kendall (Sen+M-K) trend test in this study to examine the temporal and spatial variation trends of PM_2.5 and O₃ in the Yangtze River Delta (YRD), from 2003 to 2020. We identified the regions where these pollutants exhibited synergistic changes and established the pathways between the pollutants and their potential drivers, using geographically weighted random forest algorithms and structural equation modeling. The study results revealed as follows: (1) Overall, the PM_2.5 concentrations show a decreasing trend, while the O₃ concentrations exhibit an increasing trend, in the YRD. Analysis of the combined trends indicates that approximately 95% of the area displays opposing trends for PM_2.5 and O₃, with only about 4% in the southern region showing synergistic trends for both pollutants. (2) Drought and the average temperature are the main drivers of the changes in PM_2.5 and O₃ concentrations in areas experiencing synergistic changes. Their combined effects alleviate the aggregation of PM_2.5 and reduce the formation of VOCs, indirectly reducing the generation of pollutants. The negative effect of the average temperature on the O₃ concentration may indicate the existence of nonlinear effects and complex interaction effects between the drivers. NO_x and VOCs play important dual roles in the generation and conversion of pollutants, although their overall impact is smaller than meteorological factors. They produce significant indirect effects through their interaction with meteorological and other human factors, further affecting the concentrations of PM_2.5 and O₃. In areas without coordinated changes, the main impact of meteorological factors remains unchanged, and the relationship between the two anthropogenic emission sources and their effects on PM_2.5 and O₃ are complex, with different directions and levels involved. This study provides detailed insights into the drivers of air quality changes in the YRD and offers a scientific basis for environmental management authorities to develop more comprehensive and targeted strategies for balancing the control of PM_2.5 and O₃ pollution.

1. Introduction

Air pollution events have become increasingly frequent and severe, in tandem with the rapid development of human society in recent decades. Cross-regional PM_2.5 pollution has had a significant negative impact on air quality globally [1,2]. Various regional governments in China implemented a range of measures to reduce PM_2.5 and other pollutants from 2013 to 2017, achieving some success [3]. However, PM_2.5 levels in certain areas still exceed the guidelines set by the World Health Organization (WHO) [4]. While regulatory agencies have made progress in controlling PM_2.5, O₃ pollution has become a new air quality challenge in many regions [5], especially in eastern China, where O₃ levels have risen by 1–3 ppb per year [6]. In this context, it is crucial to investigate the factors driving the synergistic changes between PM_2.5 and O₃ pollution [3,6].

Unfortunately, the research on the synergistic changes in terms of PM_2.5 and O₃ pollution and the environmental drivers of such changes is still limited, especially in regard to the large-scale spatial perspective. Most research currently focuses on the analysis of the causes of a single pollutant, that is, the study of individual pollutants such as PM_2.5 or O₃ [7,8]. Such analysis ignores the information on the interactions between the two pollutants [1] and, even if the same influencing factors are obtained, the consistent effect of these factors on the changes in terms of the two pollutants cannot be confirmed. For example, Yang et al. [9] employed the same method to separately identify the key meteorological factors affecting PM_2.5 and O₃ concentrations, but did not identify common drivers behind the variations in both pollutants. Although Wang et al. [2] established various combinations of PM_2.5 and O₃ thresholds, they did not further analyze the relationship between the two pollutants and meteorological factors across these combinations. While Wang and his colleagues utilized the correlation coefficients between the two pollutants and their relationships with meteorological factors to characterize their coordinated control, this approach did not elucidate the effects that led to the differing degrees of variation in the two pollutants.

Previous research on the synergistic variations in PM_2.5 and O₃ concentrations has predominantly relied on pollutant records from ground monitoring stations [1,2,3,5]. However, the absence of monitoring stations in certain areas or equipment malfunctions can lead to a lack of pollutant data, which in turn impedes large-scale and long-term studies [4]. Remote sensing satellite data effectively address the spatial variability in pollutant measurements [1]. Additionally, the factors influencing pollutants are quite complex, with human activities and the climate being widely recognized as key determinants of PM_2.5 and O₃ levels [2,5,10,11]. Air pollution and the climate are interlinked through a bidirectional feedback mechanism: on one hand, the climate affects the dispersion of pollutants in the atmosphere; on the other hand, air pollution can influence the balance in terms of climate development [12,13]. Meteorological conditions are external factors that significantly influence the dispersion, transportation, and deposition of pollutants [4]. For instance, temperature and precipitation are well-recognized for their impact on PM_2.5 and O₃, in this context [14,15]. In addition, several studies have shown that the self-cleaning function of green vegetation can effectively enhance air quality [4,16]. Meanwhile, biomass burning is also a significant contributor to air pollution. In the past thirty years, anthropogenic biomass burning activities have surged dramatically, contributing to the formation of tropospheric O₃, photochemical oxidants, and atmospheric brown clouds [17,18].

Given that analyzing the synergistic variations in PM_2.5 and O₃ concentrations and understanding the underlying causes of these variations has become a critical task in regard to air quality management in China, this study employed trend analysis to identify different combinations of PM_2.5 and O₃ variation degrees (hereafter referred to as PM_2.5-O₃). We then applied a combination of geographically weighted random forest algorithms and structural equation modeling to investigate the potential environmental factors influencing the synergistic (or non-synergistic) variations in PM_2.5 and O₃ concentrations. The objectives of this research are: (1) to elucidate the spatiotemporal trends and spatial distribution patterns of PM_2.5 and O₃ in the Yangtze River Delta (YRD) region between 2003 and 2020; (2) to determine the presence of regions within the YRD where PM_2.5 and O₃ exhibit synergistic variations; and (3) to explore the relationships between the pathways among the primary drivers of PM_2.5 and O₃ synergistic variation regions and their overall effects on the two pollutants. We believe that this study will significantly enhance the scientific rigor and targeted approach to air quality management in the YRD. By analyzing the synergistic variations in PM_2.5 and O₃ concentrations and their drivers, we can offer more precise policy recommendations and optimize management strategies for better air quality. Additionally, these findings provide targeted guidance for cross-regional environmental cooperation, promoting effective coordination in regard to pollution control efforts across the region.

2. Materials and Methods

2.1. Study Area

The 26 cities in the Yangtze River Delta (YRD) region of eastern China (115°54′~123°10′ E and 28°0′~34°18′ N, Figure 1a) form one of the most economically developed urban clusters in China, contributing approximately 20% of the national gross domestic product (GDP) and housing about 16% of the national population, as of 2021 [19]. The high level of human social and economic activities has long posed a severe air pollution threat to this region. According to data from the China National Environmental Monitoring Center (http://www.cnemc.cn/ (accessed on 9 May 2024)), the average concentrations of PM_2.5 and O₃ over the past decade in this region were approximately 41.2 μg/m³ (ranging from about 61.2 μg/m³ in 2014 to around 30.6 μg/m³ in 2023) and 97 μg/m³ (ranging from about 90.4 μg/m³ in 2014 to approximately 100.2 μg/m³ in 2023), respectively (Figure 1b).

2.2. Data Source and Processing

2.2.1. Air Pollutant Data (PM_2.5 and O₃)

The air pollution data used in the study comprises annual average PM_2.5 and O₃ levels from 2003 to 2020. The PM_2.5 data were provided by the Atmospheric Composition Analysis Group at the University of Washington “https://sites.wustl.edu/acag/datasets/ (accessed on 9 May 2024)” [20,21]. Because of the high spatial resolution of the dataset (0.01° × 0.01°), it has been widely used to study the health effects of environmental pollution and exposure to air pollution [4,9,22]. Van Donkelaar et al. [20] and Hammer et al. [21] confirmed the validity of these data in China, as the data exhibit a linear correlation with ground-based observations from the region. Additionally, we used total column ozone data obtained at 1000 hPa (daytime only/upward flow) from AIRS3STM version 7.0, provided by NASA. These data, which have a spatial resolution of 1° × 1°, were derived from measurements taken by the Atmospheric Infrared Sounder (AIRS) on the Aqua satellite [23]. AIRS, in combination with the Advanced Microwave Sounding Unit, constitutes an advanced atmospheric observation system that delivers enhanced measurement accuracy (Aumann et al., 2003) [24]. Unlike other sensors, this combination offers superior cloud-clearing capability (around 80%) [23,25]. The data are now extensively used in research on air pollution and health [23,24,25,26].

2.2.2. Humanistic and Social Data

Human societal development is a significant factor contributing to increased air pollution. Societal indicators, such as population density (POP), gross domestic product (GDP), the degree of human disturbance (DH), and road density (roadD), were taken into consideration in this study. POP was obtained from https://hub.worldpop.org/ (accessed on 9 May 2024). WorldPop is a global population data assessment project by the University of Southampton and its dataset includes numerous socio-economic attributes, such as population density [5]. These data, with a spatial resolution of 1 km, were estimated using a top-down random forest algorithm, based on the decomposition of global census data [27,28]. Over the years, these data have been widely applied in various research fields within Earth Sciences [5,29].

GDP data for 2000, 2005, 2010, 2015, and 2019 were provided by the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences “https://www.resdc.cn (accessed on 9 May 2024)”. These data, available as 1 km spatial grid data, were obtained from county-level GDP statistics and incorporate spatial interactions with land-use types, nighttime light intensity, and settlement density [30]. We calculated the annual GDP growth rate from the statistical yearbook and used the grid calculator in ArcGIS 10.8 to supplement the GDP for missing years [4].

Human activities and the disturbance caused create diverse ecological and environmental effects. This study quantified the intensity of human activity by calculating the Human Influence Index through a weighted overlay of land-use data, based on state changes [31]. This approach not only reduces collinearity between different land-use types [4], but also improves the understanding of the driving mechanism of human activity intensity in terms of environmental pollution changes. This insight helps to effectively regulate human activities and prevent or mitigate potential environmental crises [31]. We obtained land cover raster data, with a 300 m spatial resolution, for the time between 2003 and 2020, from the Copernicus Climate Change Service data platform “https://cds.climate.copernicus.eu/ (accessed on 9 May 2024)”. Following the weight allocation criteria proposed by Beyhan et al. [32], we assigned weights to various land-use types to quantify the degree of human disturbance on the ecosystem.

Road construction and transportation planning are two major factors influencing human activity areas, with profound effects on urban air quality [29]. Road density serves as a valuable indicator for assessing urban traffic conditions, regional socio-economic development levels, and the extent of human activities [33]. We obtained road vector data, at the 1:250,000 scale, from the National Geographic Information Directory Service “https://www.webmap.cn (accessed on 9 May 2024)”. These data were rasterized using the Line Density tool in ArcGIS 10.8, resulting in a road density raster with a 1 km resolution [16].

In addition, the grid data for NO_x, VOCs, and SO₂ were provided by Tsinghua University’s Multi-resolution Emission Inventory for China (MEIC). The data cover a number of emission sources, such as industry, energy, and transportation, and have monthly temporal resolutions and a 0.25° spatial resolution. The MEIC data are widely used in many climate and air quality studies, such as those evaluating the effect of emission reduction policies and simulating regional air quality changes [34].

2.2.3. Climatic Data

Compared to many other meteorological factors, temperature is most well-known for its impact on PM_2.5 and O₃ concentrations, by interfering with atmospheric stability [35,36]. The diurnal temperature range plays a significant role in driving changes in air pollutants and affecting air visibility, by stimulating urban heat islands [37,38]. Also, some researchers have found that drought can potentially suppress air pollutants, by affecting the physiological processes of vegetation [39,40]. Our previous studies on tropical and subtropical regions have confirmed the significant influence of these three variables on the distribution and variation in regional air pollutants [4,16]. Considering the complex correlation between air pollutant levels and the climatic factors described above, both spatially and temporally [4,40,41], we used Climatic Research Unit gridded time series (CRU TS) data, produced by the National Centre for Atmospheric Sciences (NCAS) in the United Kingdom, recognized as having high use value [29,42]. These data include grid data on the monthly average temperature, monthly average diurnal temperature ranges, monthly average precipitation, and monthly average potential evapotranspiration, during 2003 to 2020 (with a spatial resolution of 0.5°), which can be obtained from https://crudata.uea.ac.uk/cru/data/hrg/ (accessed on 9 May 2024). To maintain consistency in regard to the time scale, we used the Cell Statistics tool in ArcGIS 10.8 to synthesize the appellate data into annual meteorological parameters, including the annual mean temperature (TMP), annual mean diurnal temperature range (DTR), annual cumulative precipitation (PREC), and total annual evapotranspiration (PET). We then used the Raster Calculator to calculate the ratio of the annual cumulative precipitation (PREC) to the total annual evapotranspiration (PET) as a measure of drought, in combination with the aridity index (AI) [41]. The wind speed data came from ERA5 reanalysis data provided by the Copernicus Climate Data Store (CDS) “https://cds.climate.copernicus.eu/ (accessed on 9 September 2024)”. The data cover a global scale and provide an hourly wind speed component, with a spatial resolution of 0.25°. These high-resolution data are widely used in the field meteorology, environmental science, and engineering.

2.2.4. Vegetation Data

We selected two vegetation indices that represent vegetation growth to characterize the vitality and health of the ecosystems in the study area: one is the Normalized Vegetation Index (NDVI), which measures vegetation health and density [43]. The other is the Global Vegetation Moisture Index (GVMI) [44], which monitors vegetation health and growth, and minimizes the perturbation effects of geophysical and atmospheric influences compared to other vegetation indices [45]. In this study, we obtained the NDVI data, over the period from 2003 to 2020 (with a spatial resolution of 250 m), from the 16-day synthetic product MOD13Q1, provided by NASA. We also used the near-infrared band (NIR, 700~1000 nm, with a spatial resolution of 250 m) and short-wave infrared band (SWIR, 1000~2500 nm, with a spatial resolution of 1000 m) from MOD09A1, an 8-day synthetic product provided by NASA, to calculate the GVMI. The formula is as follows:

$$\mathrm{G}\mathrm{V}\mathrm{M}\mathrm{I}={\displaystyle \frac{\left(NIR+0.1\right)-\left(SWIR+0.02\right)}{\left(NIR+0.1\right)+\left(SWIR+0.02\right)}}$$

ENVI software (version 5.3) was used to carry out the processing of the above-mentioned data that can be obtained from https://search.earthdata.nasa.gov/ (accessed on 9 September 2024).

2.2.5. Biomass Burning Emissions Data

The global fire emissions database (GFED) “http://www.globalfiredata.org (accessed on 9 May 2024)” provides ideal data on biomass combustion emissions. At present, the data have been updated to version 4.1 and the annual emissions are divided into grassland, grassland and shrub fires, forest fires (northern, temperate, and deforestation), peat fires, and agricultural waste incineration [46]. We obtained the data on carbon and dry matter emissions from the database and the spatial resolution of the data was 0.25°. These data have been widely used in scientific research on the earth’s ecological environment [29,47,48].

2.2.6. Sample Form and Data Extraction

Considering that the spatial resolution of the data used in this study is inconsistent and that the resolution of the data on the meteorological factors and biomass burning emissions is too coarse, a very small grid scale may result in most sample units exhibiting the same attribute characteristics, which cannot effectively reflect the spatial variation of these factors. Therefore, we used the “CreateFishnet” tool in ArcGIS 10.8 to resample all the variables into a 5 km × 5 km grid. A total of 11,110 sample units were created in the YRD. Subsequently, we extracted variable data with these grids. First, we categorized the trend changes in PM_2.5 and O₃ concentrations based on the Sen+M-K method (described in detail in Section 2.3). We then overlaid the different trends in PM_2.5 and O₃ concentrations in ArcGIS to obtain five categories of PM_2.5 and O₃ trend groups (abbreviated as PM_2.5-O₃) (Figure 2). Next, we extracted the PM_2.5-O₃ type for each cell as the dependent variable. Finally, we extracted the mean values for the years 2003–2020 and the annual data for all the data from Section 2.2.2, Section 2.2.3, Section 2.2.4 and Section 2.2.5, which were used as explanatory variables for the GWRFC and SEM analyses. Table S1 and Figure S1 present the basic statistics and spatial distribution of the annual mean values of these factors, respectively.

2.3. Theil–Sen Median Trend Analysis and Mann–Kendall Test (Sen+M-K Test)

We employed Theil–Sen median trend analysis to assess the long-term trend amplitude of PM_2.5 and O₃ in the study area. It is a nonparametric statistical method commonly used in trend estimation involving long time-series data [49,50]. At the same time, the Mann–Kendall test was performed on the 18-year air pollution data collected, to test the direction of change in pollutant concentrations. We used the Theil–Sen median method to analyze the pixel-by-pixel trends in PM_2.5 and O₃ concentrations during the time span from 2003 to 2020. The specific formula is as follows [50,51]:

$$S_{pollution}=\mathrm{M}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}\left({\displaystyle \frac{{pollution}_j-{pollution}_i}{j-i}}\right),2003\le i\le j\le 2020$$

where S_pollution is the changing trend of PM_2.5 or O₃; pollution_j and pollution_i are the observed values of PM_2.5 or O₃ in the jth year and ith year, respectively; and i < j. S_pollution > 0 (or <0) indicates that PM_2.5 or O₃ shows an upward (or downward) trend in time series data.

The Mann–Kendall (M-K) test is a nonparametric statistical method, widely used in recent years to detect the trends in data time series [48]. The basic idea is to examine trends in the sequence by comparing the size relationship between each data point and the preceding data points, and then determine the direction of the trend based on the positivity or negativity of the rank sum. This method is a nonparametric method, with the advantages of being nonlinear, a non-assumption being applied in terms of normal distributed data, and non-affection being applied in regard to missing values and abnormal values [49]. The calculation process for the M-K test consists of the following steps [50,52,53]:

(1)

Define the test statistic S:

$$S=\sum _{j=1}^{n-1}\sum _{i=j+1}^{n}sgn\left({pollution}_j-{pollution}_i\right)$$

$$sgn\left({pollution}_j-{pollution}_i\right)=\left\{\begin{array}{c}1, { pollution}_j-{pollution}_i>0\\ 0, {pollution}_j-{pollution}_i=0\\ -1,{ pollution}_j-{pollution}_i<0\end{array}\right.$$

where n is the time series length of PM_2.5 or O₃ and sgn () is a symbolic function.

The statistic S roughly follows a normal distribution, and without considering the existence of equivalent data points in the sequence, its expectation and variance are [50]:

$$\mathrm{E}\left(S\right)=0, \mathrm{V}\mathrm{a}\mathrm{r}\left(S\right)={\displaystyle \frac{n(n-1)(2n+5)}{18}}$$

(2)

The statistic Z in a normalized M-K test is calculated as follows [49]:

$$Z=\left\{\begin{array}{c}{\displaystyle \frac{S-1}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}\left(S\right)}}}, S>0\\ 0, S=0\\ {\displaystyle \frac{S+1}{\sqrt{\mathrm{V}\mathrm{a}\mathrm{r}\left(S\right)}}}, S<0\end{array}\right.$$

At a given significance level α = 0.05, when |Z| > 1.96, the PM_2.5 or O₃ time series data have a significant trend, while the trend was not significant on the contrary.

Combined with the Theil–Sen median trend analysis and Mann–Kendall test results, |Z| > 1.96 passed the significance test with a confidence level of 95%. We classified the interannual variation trends in PM_2.5 and O₃ concentrations into five categories: significant deterioration (S < −0.0005, Z < 1.96); mild deterioration (S < −0.0005, −1.96 ≤ Z ≤ 1.96); essentially unchanged (−0.0005 ≤ S < 0.0005, −1.96 ≤ Z ≤ 1.96); mild improvement (S ≥ 0.0005, −1.96 ≤ Z ≤ 1.96); and significant improvement (S ≥ 0.0005, Z > 1.96) (Liu et al., 2020) [54].

We used the “terra”, “trend”, and “dplyr” packages in the R environment to visualize the results of the Sen+M-K test.

2.4. Geographically Weighted Random Forest Classification (GWRFC)

The geographically weighted random forest (GWRF) is a machine learning algorithm used to analyze and explore spatial data, which combines the advantages of both the GWRF and random forests (RFs). The GWRF can explain the spatial variation relationships between the response variable and the predictive factors, while accounting for the nonlinear and interaction effects of the independent variables [55]. This algorithm first divides the study area into a set of smaller regions, with separate models being constructed using an RF for each region. Then, it estimates the spatial variation relationships between the variables and the independent variables within each area [56]. In our research, the parameters of the GWRF are set according to the RF. We performed the calculations using the “GWRFC” package in the R environment, which uses the “ranger” package for parallel computation of the RF to improve the processing speed [57]. The GWRFC can generate the importance of the variable, classification probabilities, and accuracy of RF models at the local level.

We used the kappa coefficient both for consistency testing and for measuring the explanatory ability of the variables used in the modeling method for PM_2.5-O₃. The kappa coefficient is a testing metric for multiclass data, based on the confusion matrix. It is used to measure classification accuracy and is primarily used for consistency testing [58]:

$$kappa={\displaystyle \frac{p_0-p_e}{1-p_e}}$$

$$p_0={\displaystyle \frac{{\sum }_{i=1}^n{a}_{ii}}{N}}$$

$$p_e={\displaystyle \frac{{\sum }_{i=1}^n{a}_{i+}{b}_{i+}}{N^2}}$$

where p₀ represents the proportion of observed accuracy or consistency, which indicates the overall classification accuracy. Moreover, p_e denotes the proportion of units that are expected to be consistent by chance or the proportion of random consistency. Where a_ii represents the elements on the diagonal of the confusion matrix. Then, a_i+ and b_i+ denote the sum of all the elements in the ith row and ith column, respectively. N represents the total number of samples. The kappa value typically ranges from 0 to 1. In this study, it was categorized into six ranges to represent different levels of consistency: 0.0~0.20 for slight consistency, 0.21~0.40 for fair consistency, 0.41~0.60 for moderate consistency, 0.61~0.80 for substantial consistency, 0.81~0.90 for high consistency, and 0.91~1 for almost perfect consistency.

2.5. Structural Equation Model (SEM)

Considering that we need to analyze PM_2.5 and O₃ as dependent variables simultaneously, to understand how they are influenced by human and climate factors, we used the “lavaan” and “semPlot” packages in R studio to construct the structural equation model (SEM) and generate path diagrams [59,60]. The development of an SEM is a multivariate analysis method used to validate the relationships between one or more independent variables and one or more dependent variables [61]. The core idea is to analyze the relationships between the variables based on their covariance matrix, serving as a confirmatory statistical method. Essentially, it is a generalized form of linear modeling [61,62]. Structural equation modeling integrates two statistical techniques, namely factor analysis and path analysis. It allows for the simultaneous consideration and handling of multiple dependent variables and their relationships with predictive variables and helps in constructing a causal network of these dependent variables to build a framework model [29,61,63,64]. Structural equation modeling has recently been widely used in environmental ecology [29,59,60]. Based on previous research, the model’s validity was assessed using the following goodness-of-fit indices: the root mean square error of approximation (RMSEA < 0.05), standardized root mean square residual (SRMR < 0.08), goodness-of-fit index (GFI > 0.90), comparative fit index (CFI > 0.90), and non-normed fit index (NNFI > 0.90) [59,63].

3. Results

3.1. Change Trend Analysis of PM_2.5 and O₃

By overlaying the Theil–Sen median trend classification results with the Mann–Kendall test classification results, we obtained the distribution maps of the PM_2.5 and O₃ concentration trends (Figure 2a,b). The figure depicts the overall trends and significance of the PM_2.5 and O₃ concentration changes in the YRD throughout the years 2003 to 2020. The overall trend in the PM_2.5 concentration in the YRD showed approximately 73.3% of the region experiencing a significant improvement (Figure 2a). Only localized areas in the eastern and northeastern parts showed insignificant and clustered trends. This indicates that the PM_2.5 control measures in the study area over the past 20 years have been effective. Unlike PM_2.5, the O₃ concentrations showed almost a rising trend uniformly, with little difference between the two worsening categories (MD: 51.69%; SD: 45.64%, Figure 2b). In addition, after overlaying the change trends of the two pollutants, we identified five combinations of PM_2.5-O₃ change relationships (Figure 2c). The combination where PM_2.5 significantly decreased and O₃ slightly increased (abbreviated as PM_2.5(−2)-O₃(1)) occupied the largest proportion, nearly half (46.43%) of the study area. These regions were primarily located in the southern part of the area. Moreover, 24.15% of the study area experienced a significant increase in O₃, while PM_2.5 significantly decreased (abbreviated as PM_2.5(−2)-O₃(2)), primarily concentrated in the central and northern regions. This was followed by regions (approximately 20.81% of the study area) where PM_2.5 slightly decreased and O₃ significantly increased (abbreviated as PM_2.5(−1)-O₃(2)), which were mainly located in the eastern and northeastern parts of the area. The area with a slight improvement in PM_2.5 and a slight increase in O₃ (abbreviated as PM_2.5(−1)-O₃(1)) accounted for only 5.15% of the study area. Moreover, only a small proportion (about 3.5%) of the region experienced a decreasing trend in both pollutants (abbreviated as PM_2.5(−1, −2)-O₃(−1)).

3.2. Performance of the GERFC

We determined the optimal bandwidth for the model, based on the kappa coefficient (Figure S2). As the bandwidth rose from 20 to 260, the mean kappa index rose from 0.859 to 0.961 and the standard deviation decreased. However, within the bandwidth range of 260–300, both the mean kappa coefficient (ranging from 0.958 to 0.961) and the standard deviation remained nearly constant (0.042–0.047). Therefore, we conducted GWRFC fitting for five types of PM_2.5-O₃ change patterns and all the independent variables within the bandwidth limit of 260. As shown in the predicted probability distribution diagram for the five categories of PM_2.5-O₃ (Figure S3), all the categories exhibited high predicted probabilities (>0.95). Figure 3 illustrates the consistency of the complete variable samples and classified variable samples between the predicted categories (as calculated by the GWRFC with a bandwidth of 260) and the actual observed categories. Whether using the four-category variable datasets or the complete variable datasets, the kappa coefficient for all the samples almost exceeded 0.6 (Figure 3a–e). This suggests that the latent factors, discussed in Section 2.2, have substantial explanatory power in regard to PM_2.5-O₃, particularly the artificial and meteorological variables (Figure 3b) and the complete variable set (Figure 3e), where the results are nearly identical. The worst consistencies were seen in regard to fire emissions (Figure 3d) and vegetation factors (Figure 3c). From Figure 3f, we can observe that approximately 91.1%, 91.5%, and 86.4% of the samples in the complete, artificial, and meteorological variable datasets, respectively, fell within the range of 0.9 to 1. In contrast, about 45.8% and 43.2% of the samples in the vegetation and fire emission variable datasets, respectively, were within this range. Obviously, human factors exhibit the greatest explanatory power and broadest scope in terms of PM_2.5-O₃ compared to the other three types of variables, with meteorological factors coming next.

Additionally, we divided the complete dataset into five subsamples, fitted each subsample separately, and assessed the robustness of the GWRF based on the model’s kappa index and the ranking of the variables’ importance. As shown in Figure S4, the mean, median, and quartile ranges of the kappa index of the five subsamples were almost the same, falling within the range of 0.9 to 1. On the other hand, we also obtained the variables’ importance results for each subsample. The variables that had a significant impact on PM_2.5-O₃ were mainly concentrated in the three meteorological factors of AI, DTR, and TMP, as well as the two human factors of GDP and DH (Figure S5).

3.3. The Importance of the Influence Factor for PM_2.5-O₃

Based on the results of the local variable importance, we presented the statistical measures of the variable importance scores (Table S2) and displayed the spatial distribution of all the variable importance scores (Figure 4). AI, DTR, and TMP have higher importance scores and are more widely distributed (Figure 4a–c). Their statistics were obviously higher than the statistics for the other variables (Table S2). AI showed lower scores in parts of the northeast and south, while other areas exhibited higher scores. The distribution pattern for DTR is similar, but its score is slightly lower than AI. High scores for TMP were predominantly concentrated in the southern, western, northwest, and northern regions. GDP and DH are two crucial anthropogenic factors influencing PM_2.5-O₃ variations. The distribution of GDP and DH scores is similar, but there are still differences. The distribution of GDP leans toward the central and western regions, while DH is more evenly distributed. The importance-ranking chart for the explanatory variables for different types of PM_2.5-O₃ changes further confirmed the significance of the five variables mentioned above (Figure 5). We also found that TMP, DTR, and AI were all important explanatory variables for any type of PM_2.5-O₃ change and the importance scores for the three factors were relatively similar.

To understand the changes in the variable importance ranking in different spatial regions intuitively, we mapped out the spatial distribution of the variables’ importance ranking (Figure S6). We found that in regard to the spatial distribution of the top five variable importance rankings, human factors and meteorological factors accounted for more than 65% of the total, which were 83.8%, 85.3%, 80.8%, 74.4%, and 67.4% in turn. In the top ranking, TMP, DTR, and AI account for a relatively high proportion of the area; then, the influence of human factors begin to stand out, especially GDP. The importance of three artificial emission factors also begin to stand out.

3.4. Path Analysis of Major Influencing Factors on PM_2.5 and O₃

Based on the results detailed in Section 3.1, we divided the study area into regions with synergistic changes in PM_2.5 and O₃ concentrations and regions without such changes. We then performed path analysis on the annual average PM_2.5 and O₃ data for both types of regions, as well as on the top five most important influencing factors identified in Section 3.3 (i.e., the factors selected by the GWRFC). To achieve an appropriate model fit, we adjusted the paths between the variables. The resulting path distribution map is plotted in Figure 6, with the model fitting indicators all falling within the specified range (for regions with synergistic changes: GIF = 0.996, CIF = 1, NNIF = 1.001, RMSEA = 0.000, SRMR = 0.012; for regions without synergistic changes: GIF = 1, CIF = 1, NNIF = 0.998, RMSEA = 0.027, SRMR = 0.001).

In the region with PM_2.5 and O₃ co-changes, the explanatory variables in regard to the pathways had strong explanatory power for PM_2.5 and O₃ changes (R² is 0.89 and 0.87, respectively) (Figure 6a) and all the pathways were significant (p-value < 0.05). We found that AI had the largest direct negative impact on PM_2.5 (normalized coefficient of −0.861), followed by DTR and VOC, with direct positive effects of 0.174 and 0.169, respectively. Although TMP had a small direct effect on PM_2.5, it still maintained a significant negative correlation (−0.056). For O₃, AI had the largest direct impact (−0.757), followed by TMP, which had a direct impact of −0.315 on O₃, showing a significant temperature control effect. GDP and DTR also have a significant direct negative impact on O₃. Although the direct effect of DH on PM_2.5 is significant, its direct effect on O₃ is relatively weak. In this region, the correlation coefficient between PM_2.5 and O₃ was 0.488, indicating that the positive correlation between the two was significant. In the non-covariant region, the explanatory variables in regard to the pathway also have high explanatory power for PM_2.5 and O₃ (R² is 0.90 and 0.89, respectively) (Figure S7a) and all the factors in the pathway are significant (p-value < 0.05). DTR has the most prominent direct positive effect on PM_2.5 (0.444), while AI has a significant direct negative effect on PM_2.5 (−0.429). The direct influence of DH on PM_2.5 and O₃ in this area is remarkable and the standardized coefficients are 0.286 and −0.23, respectively. Compared with the co-varying area, AI has less direct influence on PM_2.5 and O₃. For O₃, AI is still the most important direct factor of influence, followed by DH and meteorological conditions (DTR and TMP). Different from the synergistic region, NO_x and VOCs have relatively little influence on O₃, showing the regional differences in terms of influence factors. The relationship between PM_2.5 and O₃ is slightly weak in this region (the correlation coefficient is 0.363), but it still shows a positive correlation trend.

In addition to the direct effects of the drivers on PM_2.5 and O₃ mentioned above, we also explored the indirect effects of these drivers on the two pollutants from the pathway diagrams and calculated their total effects on the two pollutants (Figure 6b,c and Figure S7b,c). In the region of synergistic change, the analysis of the indirect effects reveals some interesting conclusions. DH had the largest indirect effect on PM_2.5, although it had no direct effect on PM_2.5, suggesting other pathway effects of human activity. The indirect effect of TMP on PM_2.5 was −0.474, indicating that it played an important role in the temperature change. The indirect effect of AI on PM_2.5 was slightly lower than that of TMP (−0.425), but its total effect was still the largest among all the factors. There was no direct effect on NO_x, but the indirect effect was 0.195. The indirect effect of VOCs on PM_2.5 is very small, mainly through direct effects. The indirect effects of both DTR and GDP are weak. Among the mechanisms of O₃ influence, AI exhibited the most significant indirect and total effects, highlighting the important role of dryness in reducing O₃ concentrations. The indirect effect of TMP was −0.422 and the total effect was −0.57, indicating that air temperature had an inhibitory effect on O₃. Although the indirect effect of DH on O₃ was 0.327, the total effect reduced to 0.216, due to the direct negative effect. At the same time, the indirect effects of NO_x, DTR, and VOCs were similar (−0.19, 0.167, and −0.115, respectively). These results illustrate the complex interactions between these factors and O₃. In the non-synergistic change region, the indirect and total effects of AI on PM_2.5 were the highest, which were −1.208 and −0.78, respectively. The indirect effects of TMP, DTR, and DH were −0.459, 0.39, and 0.365, and the direct effects were −0.32, 0.834, and 0.651, indicating that the effects of DTR and DH on PM_2.5 exceeded that of TMP, NO_x, and VOCs. The indirect effects of NO_x and VOCs on PM_2.5 were 0.13 and 0.014, respectively, and the total effect of NO_x was greater than that of VOCs. For O₃, the indirect effect of DH was the maximum at 0.881, followed by DTR (0.749) and TMP (−0.401). The indirect effect of AI is smaller (−0.163), but its total effect is still the largest. The indirect effects of NO_x and VOCs on O₃ were weak, which were much lower than those in the synergistic change area.

4. Discussion

4.1. A Broad Trend of Decreasing PM_2.5 and Increasing O₃

We found the spatial distribution trend of PM_2.5 reduction in the YRD between 2003 and 2020, which is consistent with previous studies [3,5]. We observed a slight decreasing trend in the eastern, northeastern, and northern regions. One possible reason is that these areas are home to several major economic and industrial centers, including Shanghai, Changzhou, Wuxi, Suzhou, and Yangzhou (Figure 1a), which have high population densities (Figure S1d) and per capita GDP distributions (Figure S1e). Additionally, the level of human activity interference is significantly higher than in other regions (Figure S1f). Our previous research in the YRD also supports this conclusion [4]. Different from PM_2.5, we noticed the increasing trend in terms of O₃ concentrations in the study area, which is the same as the research results in other studies [3,5,6], and the increasing trend in O₃ concentration was significant in areas where PM_2.5 was slightly reduced, even in the west, south, and southeast, where the PM_2.5 concentration was significantly reduced.

We categorized five combinations of PM_2.5-O₃ trends and found that the synergistic reduction of the two pollutants only occurred in a small area in the southern part of the YRD (Figure 2c). In over 95% of the study area, there was no synergistic change between the two pollutants; specifically, the PM_2.5 concentration decreased, while the O₃ concentration increased. Data from the monitoring stations further confirm this trend. (Figure 1b). This may indicate that the pollutants’ changes are not synergistic, despite the same degree of change occurring. The reduction in PM_2.5 may be related to the reduction in local industrial emissions and transportation emission reduction policies [65] or other environmental protection measures [66]. The significant growth in O₃ may come from a combination of environmental conditions, including changes in NO_x and volatile organic compound (VOC) emissions [10], as well as changes in environmental conditions, such as increased ultraviolet radiation intensity [67], the atmospheric oxidation capacity, or oxidant concentration [68]. A significant decrease in PM_2.5 and a slight growth in O₃ do not represent a synergistic change. The significant decrease in PM_2.5 may be due to the same reasons mentioned above, while the slight increase in O₃ may be mainly due to indirect effects, such as the complexity of chemical reactions [68] and the effects of photochemical reactions [10]. For example, increased light intensity or changes in atmospheric chemicals may affect the rate of O₃ formation [67,68]. Similarly, a slight decrease in PM_2.5 and a significant increase in O₃ also indicate that the changes in the two pollutants are uncoordinated. The slight reduction in PM_2.5 may be regarded as a positive environmental improvement effect [3,5], but the significant increase in O₃ reflects that it may be necessary to pay more attention to the emission control of its precursors and the comprehensive impact of changes in atmospheric environmental conditions. In the area where PM_2.5 decreased slightly and O₃ rose slightly, despite a slight change in both pollutants, the reason for it may not necessarily explain the synergistic change between the two pollutants. Slight changes may be more of a reflection of the interfering effects between different substances in relation to environmental factors [10,69]. Studies have found that even if PM_2.5, NO_x, and VOC concentrations do not change significantly, O₃ concentrations may still vary slightly due to changes in light, temperature, and other environmental conditions [10]. The areas where both PM_2.5 and O₃ decreased slightly may be the only type of co-variation observed in this study, but such areas account for less than 4% of the total area. Since the beginning of the 21st century, China has been paying more attention to PM_2.5 emissions than to O₃ pollution. Additionally, increased vehicle exhaust fumes in cities have led to rising O₃ levels. Therefore, it is understandable that there is a widespread trend in the YRD, where the PM_2.5 concentration decreases, while the O₃ concentration increases [6].

4.2. The Drivers of Change Trends in PM_2.5-O₃

The synergistic variation of the two pollutants occurred mainly in the southwest of Hangzhou and Jinhua (Figure 2c). The PM_2.5 and O₃ concentrations in this specific region exhibited different pathway relationships in response to different environmental factors, which may indicate the complex interaction of environmental factors in the region. In this region, drought has the greatest impact on PM_2.5 and O₃ concentrations, further confirming the crucial role of dryness in regulating atmospheric pollutant concentrations. The significant negative effect of drought indicates that it helps to reduce the concentrations of PM_2.5 and O₃, which may be related to the enhanced atmospheric diffusion capacity in dry conditions [70]. The negative effect of the average temperature indicates its promotional effect on the diffusion of pollutants, which is consistent with the theory that a temperature rise promotes the thickening of the atmospheric boundary layer, thereby improving the diffusion capacity of PM_2.5 [71]. For O₃, the strongest effect of the average temperature may be the result of meteorological conditions and chemical reactions. It is common to promote the rate of chemical reactions at higher temperatures for O₃ formation, but this is not a determining condition [10]. Our unusual results may be explained as follows: first, the influence of the average temperature on O₃ production may have a nonlinear effect [72]. There may be an inverse effect of the average temperature in some ranges on the formation of O₃. For example, high temperatures may promote the rate of some chemical reactions, especially those involving NO_x. High temperatures may accelerate the formation of NO_x and its reaction with O₃ when enhanced by ultraviolet radiation in the stratosphere, thereby reducing O₃ concentrations [67,72,73]. Second, there may be a complex interaction between the average temperature and other factors (such as human interference, drought, etc.) [12,74], and these interactions may lead to the observed negative correlation effect (Figure 6a). Existing research also supports the importance of drought and temperature [40,75]. Demetillo and his colleagues confirmed the role of drought in reducing VOCs (O₃ and PM_2.5 precursors) [40]. Shikhovtsev and his colleagues also confirmed that many factors influence the variability in particulate matter concentrations. In addition to meteorological conditions, pollutant emission sources are also the main influencing factors and the interaction between them significantly affects the variation in particulate matter concentrations [75].

Although NO_x and VOCs have less direct effects on both pollutants than AI and TMP, they have an important indirect effect on PM_2.5 and O₃ through their interaction with other variables. NO_x is a precursor of O₃, which affects the generation and consumption of O₃ through complex chemical reactions. NO_x may also further act on O₃ concentrations and PM_2.5 by influencing other variables, such as VOCs and meteorological factors, despite its weak direct impact on both pollutants. Similarly, although the indirect effect of VOCs on PM_2.5 is small, its indirect effect on O₃ shows that it plays an important role in atmospheric chemical reactions. As an important precursor of O₃, VOCs generate O₃ through photochemical reactions and its indirect effects further reveal the influence of its interaction with other variables on pollutant concentrations. Existing studies also support the complex role of NO_x and VOCs. For example, Yao et al. [76] found that NO_x and VOCs form O₃ through photochemical reactions in the urban atmosphere and pointed out that a reduction in NO_x does not always lead to a decrease in O₃, because the ratio of NO_x to VOCs plays a key role in O₃ generation. On the other hand, NO_x may further affect the formation of O₃ by influencing the oxidation process of VOCs. The results of this study further extend this finding by quantifying the indirect effects of NO_x and VOCs using a SEM, revealing their importance in the interaction of different variables. In addition, the role of VOCs in atmospheric chemical reactions is not limited to O₃ generation, but may also indirectly affect PM_2.5 or O₃ concentrations by influencing the meteorological conditions [76]. For example, high temperature conditions may accelerate the volatilization and oxidation of VOCs, thereby increasing O₃ production, while high NO_x conditions may lead to O₃ depletion.

O₃ formation is strongly dependent on photochemical reactions with NO_x and VOCs [77]. NO_x and VOCs are also precursors of PM_2.5, contributing to the formation of nitrates and secondary organic aerosols (SOAs) [78]. In this study, structural equation modeling revealed that anthropogenic emissions of NO_x directly influence O₃, but only indirectly impact PM_2.5. This indirect correlation stems from the dependence of PM_2.5 formation from NO_x oxidation on atmospheric oxidant concentrations, including radicals and photooxidants [79]. In the YRD, major industrial activities include metal fabrication and machinery manufacturing, pharmaceutical and chemical industries, textile production, steel smelting, and coal-fired power plants. These sources exert significant influence on PM_2.5 and NO_x emissions and, thereby, on air quality. The PM_2.5 generated during metal fabrication and machinery manufacturing not only increases the particulate concentrations in the atmosphere, but also interacts with VOCs released by the pharmaceutical, chemical, and textile industries, forming secondary PM_2.5 through photochemical reactions [80]. Steel smelting and coal-fired power plants are primary sources of NO_x and their photochemical reactions with VOCs not only generate O₃, but may also influence the oxidation processes of VOCs, thereby augmenting the PM_2.5 concentration [81]. Under the impetus of various recent emission reduction measures and environmental policies, significant reductions in NO_x and PM_2.5 emissions have been achieved through the installation of high-efficiency desulfurization and denitrification equipment and the improvement of combustion technologies [34,80]. This, to some extent, explains the synergistic reduction in PM_2.5 and O₃ concentrations in response to the reduction in NO_x emissions.

In the region without synergistic change, drought has a high impact on both pollutants. The direct and indirect effects of the daily temperature range and human interference on PM_2.5 and O₃ concentrations have been shown to cause significant enhancements (Figure S7). The daily difference in average temperature may affect the pathway of O₃ generation in photochemical reactions [82], while its impact on PM_2.5 may be more due to its influence on atmospheric stability. In addition, daily temperature fluctuations may indirectly alter the PM_2.5 concentration by affecting the humidity and wind speed in the air [83]. Therefore, regions with a trend of collaborative change may exhibit more pronounced correlations in terms of the consistency of the influencing factors, while regions without collaborative change may involve more complex and diverse influencing mechanisms [2,84]. This distinction is of great significance for understanding environmental responses in different regions and formulating targeted environmental policies. Generally, a high level of human interference means more emissions from artificial activities, all of which directly increases the concentration of pollutants in the atmosphere. However, the impact of the two individual emissions on PM_2.5 and O₃ is relatively weak. This may be due to the combined effects of regional specificity, the complexity of the meteorological conditions, the secondary generation processes of the pollutants, implemented pollution control measures, and limitations in terms of the data and models.

4.3. Limitation

Combining variable selection with the GWRFC and path relationship analysis using structural equation modeling yielded excellent results in this study. Nevertheless, while the satellite data used can address the spatial coverage limitations of ground monitoring stations by providing continuous spatial information on PM_2.5 and O₃, its temporal resolution is generally lower than that of real-time data from ground stations. As a result, this study could not analyze the effects of drivers on these pollutants at a finer temporal scale. Future research will consider using daily-scale data from monitoring stations or integrating field-sampling data to address this issue. Secondly, we acknowledge potential limitations related to both the data and the model. There may be other factors that have not been taken into account, such as the temperature-disturbed atmospheric boundary layer we already mentioned, which will be a focus of our future research. On the other hand, the interaction between NO_x and VOCs may have an effect on human activities. We will further explore this complex interaction in future research to develop more effective pollution control strategies. The lack of advanced algorithms for in-depth application of multi-criteria analysis may limit our comprehensive understanding of the complex relationships and synergies between industrial pollution sources. In the next step, we will focus on applying advanced algorithms in multi-criteria analysis, such as the fuzzy comprehensive evaluation method and multi-objective linear rule, to delve into the detailed mechanisms of the interactions and photochemical reactions between different industrial pollution sources. This will help us to better quantify the impact of each source on air quality.

5. Conclusions

In this study, we used the Sen+M-K trend analysis to reveal the overall decreasing trend in PM_2.5 concentrations and the overall increasing trend in O₃ concentrations in the YRD, between 2003 and 2020. While the authorities effectively control PM_2.5 levels, O₃ concentrations have been rising, suggesting that the trends and patterns of different types of pollutants may differ. We then overlaid the trend changes of the two pollutants and discovered that in approximately 95% of the YRD region, the trends for PM_2.5 and O₃ were the opposite. This could reflect the independent effects of different pollution sources or environmental factors on these pollutants. However, in nearly 4% of the southern region, we observed that PM_2.5 and O₃ exhibited a concurrent trend, suggesting that shared drivers influence pollutant changes. We employed the GWRFC and an SEM to identify and analyze the pathways through which common drivers affect the two pollutants. Our analysis revealed that in the synergistic change region, drought has the greatest negative effect on PM_2.5 and O₃ concentrations. This suggests that dry conditions help to reduce the concentration of both pollutants at the same time, indicating a possible connection with the enhanced atmospheric dispersion capacity in dry environments. The negative effect of the average temperature on PM_2.5 concentrations further proves that a temperature increase can promote the thickening of the atmospheric boundary layer and improve the diffusion ability of pollutants. However, the response of O₃ to the mean air temperature shows an anomalous negative correlation, which may be due to the nonlinear effect of temperature on O₃ generation and the complex interaction with other potential influencing factors, which needs to be further confirmed in subsequent studies. Although the direct effects of NO_x and VOCs on PM_2.5 and O₃ are not as significant as meteorological factors, they have a significant indirect impact through their interaction with meteorological and other anthropogenic factors, further affecting the concentrations of PM_2.5 and O₃. These complex mechanisms show that NO_x and VOCs play a dual role in the process of pollutant generation and transformation. In the non-synergetic change area, meteorological factors play a leading role in the concentration of the two pollutants and human interference has a more prominent impact on O₃ and PM_2.5 concentrations. The opposite mechanism of NO_x and VOCs in regard to the two pollutants also shows the complexity of the driving factors of the two pollutants in these areas as a result of artificial emission sources. The findings of our study offer refined insights into the factors driving changes in the air quality in the YRD. This research, by exploring synergistic and non-synergistic pollutant drivers, provides a scientific basis for environmental management agencies to create more comprehensive air quality policies and balanced control strategies for PM_2.5 and O₃ pollution across various regions.