Insighting Drivers of Population Exposure to Ambient Ozone (O3) Concentrations across China Using a Spatiotemporal Causal Inference Method

Abstract

Ground-level ozone ($O_3$) is a well-known atmospheric pollutant aside from particulate matter. China as a global populous country is facing serious surface $O_3$ pollution. To detect the complex spatiotemporal transformation of the population exposure to ambient $O_3$ pollution in China from 2005 to 2019, the Bayesian multi-stage spatiotemporal evolution hierarchy model was employed. To insight the drivers of the population exposure to ambient $O_3$ pollution in China, a Bayesian spatiotemporal LASSO regression model (BST-LASSO-RM) and a spatiotemporal propensity score matching (STPSM) were firstly applied; then, a spatiotemporal causal inference method integrating the BST-LASSO-RM and STPSM was presented. The results show that the spatial pattern of the annual population-weighted ground-level $O_3$ ($PWGL{O}_3$) concentrations, representing population exposure to ambient $O_3$, in China has transformed since 2014. Most regions (72.2%) experienced a decreasing trend in $PWGL{O}_3$ pollution in the early stage, but in the late stage, most areas (79.3%) underwent an increasing trend. Some drivers on $PWGL{O}_3$ concentrations have partial spatial spillover effects. The $PWGL{O}_3$ concentrations in a region can be driven by this region’s surrounding areas’ economic factors, wind speed, and $PWGL{O}_3$ concentrations. The major drivers with six local factors in 2005–2014 changed to five local factors and one spatial adjacent factor in 2015–2019. The driving of the traffic and green factors have no spatial spillover effects. Three traffic factors showed a negative driving effect in the early stage, but only one, bus ridership per capita (BRPC), retains the negative driving effect in the late stage. The factor with the maximum driving contribution is BRPC in the early stage, but $P{M}_{2.5}$ pollution in the late stage, and the corresponding driving contribution is 17.57%. Green area per capita and urban green coverage rates have positive driving effects. The driving effects of the climate factors intensified from the early to the later stage.

1. Introduction

Ground-level ozone ($GL{O}_3$) is a well-known major atmospheric pollutant threatening human health, vegetation and biodiversity worldwide [1,2,3]. $GL{O}_3$ is also an important greenhouse gas contributing to climate change [4]. Significant increasing trends of annual $GL{O}_3$ concentrations occurred in most global urban areas from the 1990s, particularly from 2005 to 2014, and annual $GL{O}_3$ concentrations increased at a fast rate in East Asia, North America, and Europe [5]. China has established several air pollution control policies and programs [6,7] since 2016. Consequently, the concentrations of most air pollutants, including $P{M}_{2.5}$, $N{O}_2$ and $S{O}_2$, of China have decreased since 2013, but $GL{O}_3$ concentrations have shown an increasing trend [8]. $GL{O}_3$ may have become the second most hazardous pollutant after $P{M}_{2.5}$ [9]. Lyu et al. [10] pointed out that the ratio of the population living in regions with high $P{M}_{2.5}$ pollution decreased from 2015 to 2021; however, this ratio remained almost unchanged for the ozone and even increased.

Considering the complex $GL{O}_3$ concentration problems in China, a deep understanding of the spatiotemporal trends of $O_3$ concentrations in China can provide helpful evidence for environmental epidemiology and appropriate regional control strategies of $GL{O}_3$. Some researchers have investigated the spatiotemporal characteristics and influencing factors of $GL{O}_3$ in a specific city or a limited spatial region in a short-term period, mainly using the in situ monitored $O_3$ concentrations data. For instance, An et al. [11] explored the temporal variations of $O_3$ concentrations in Nanjing, Eastern China, and Wang et al. [12] studied the spatial and temporal patterns of observed ground-level $O_3$ concentrations in China. The above research found $O_3$ concentrations to be highest in summer and lowest in winter, with significant diurnal variations. Zhang et al. [13] used the spatial–temporal kriging model to analyze the daily maximum 8 h ozone concentration data from 55 air quality monitoring stations in China’s Pearl River Delta Metropolitan Region in 2015, revealing the $O_3$ concentrations have a spatial autocorrelation that undergoes temporal migration and conversion across the seasons. Tang et al. [14] and Gao et al. [15] underscored the influence of anthropogenic emissions, particularly in urban areas, on ozone levels. Liu et al. [16] suggested that weather patterns and regional interactions should be considered when designing ozone mitigation strategies. Lyu et al. [17] demonstrated the use of machine learning algorithms for predicting ozone levels, showing promising results. These findings emphasize the need for region-specific and season-specific strategies to effectively address the issue of ozone pollution in China. Some scholars have explored the space-time variations of the $O_3$ concentrations across China using simple descriptive statistical methods based on in situ monitored air quality data. Yang et al. [18], Xue et al. [19], and Li et al. [20] explored the temporal variation and spatial distribution of $GL{O}_3$ concentrations in China during 2013–2017, 2005–2017, and 2015–2020, respectively. Based on multi-source data and machine learning, Liu et al. [21] and Wei et al. [22] estimated and investigated briefly the temporal and spatial trends of the $GL{O}_3$ concentrations in China during 2005–2017 and 2013–2020. The spatial distribution of the population was not considered in the above research. However, the health exposure risk of $GL{O}_3$ concentrations is determined by $GL{O}_3$ concentrations and population spatial distribution jointly. Given this, our study focused on the population exposure to the $GL{O}_3$ concentrations, instead of semplice $GL{O}_3$ concentrations.

The formation regime or driver of $GL{O}_3$ concentrations is also a vital issue. A large number of studies have investigated the complicated relationship between $GL{O}_3$ production and main precursors, volatile organic compounds (VOC) and ${NO}_x$, from the perspective of atmospheric physicochemistry [23,24,25,26,27]. Some scholars have researched the spatiotemporal characteristics of $GL{O}_3$ formation regime of local areas of China, such as Beijing [28,29,30] and Shanghai [31,32]. Li et al. [20] and Lu et al. [33] studied the spatiotemporal variations of $GL{O}_3$ formation regimes across China. Although it is well known that meteorological factors have an influence on $O_3$ concentrations, the specific influencing mechanism still needs to be further explored. The meteorological factors’ influencing effects on $GL{O}_3$ concentrations were studied in some articles [34,35,36,37,38]. These studies all found that temperature was the dominant driver of $GL{O}_3$ concentrations at the annual time scale.

Most of the existing research has focused mainly on the pure $GL{O}_3$ concentrations without considering the population spatial distribution, and the study phases of these studies were within a short time range. Additionally, to our knowledge, few studies have investigated the simultaneous influence of natural and anthropogenic factors on population exposure to $GL{O}_3$ concentrations across China. Given this, our paper utilised the advanced Bayesian spatiotemporal hierarchy model to investigate transitions of spatial distributions of the annual $PWGL{O}_3$ concentrations over China in recent 15 years, and then presented a spatiotemporal causal inference method for the multiple-factor analysis scenario, which was used to identify the potential drivers of population exposure to $GL{O}_3$ concentrations.

2. Methods

2.1. Variable and Data

2.1.1. Population Exposure to $GL{O}_3$

The outcome variable in this study is population exposure to $GL{O}_3$ concentrations. In this study, population exposure to $GL{O}_3$ concentrations is represented by the annual population-weighted $GL{O}_3$ ($PWGL{O}_3$) concentrations calculated by the grid data of the $GL{O}_3$ concentrations and population distribution by the following formula,

(1)$PWGL{O}_{3,it}=\frac{{\sum }_{k=1}^{I_{ik}}GL{O}_{3,ikt}*Po{p}_{ikt}}{{\sum }_{k=1}^{I_{ik}}Po{p}_{ikt}}$

where $GL{O}_{3,ikt}$ and $Po{p}_{ikt}$ represent the $GL{O}_3$ concentrations and distributed population count of the $k^{th}$ grid with spatial resolution of $10\mathrm{k}\mathrm{m}\times 10\mathrm{k}\mathrm{m}$ in the $i^{th}$ prefecture region. $I_{ik}$ is the total number of the 10 km grid of across the $i^{th}$ prefecture region. The grid data of the $GL{O}_3$ concentrations at a spatial resolution of $10\mathrm{k}\mathrm{m}\times 10\mathrm{k}\mathrm{m}$ was collected from the CHAP dataset (https://weijing-rs.github.io/product.html (accessed on 6 June 2022)), which was estimated using the extended space–time extra-trees (STET) model [21] from big data including ground-based observations, remote sensing products, atmospheric reanalysis, and emission inventory [21]. The $GL{O}_3$ estimates are highly consistent with the ground-based measurements with $R^2=$ 0.96 and the root-mean-square error (RMSE) of 8.64 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$. Additionally, the $GL{O}_3$ estimates for the whole of China are highly consistent with surface measurements ($R^2$ = 0.87), and the values of the mean RMSE, the mean absolute error (MAE), and the mean relative error (MRE) over the entire domain, were 17.10 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$, 11.29 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$, and 18.38%, respectively [21]. The grid data of the population density with the spatial resolution of $1\mathrm{k}\mathrm{m}\times 1\mathrm{k}\mathrm{m}$ were collected from an open database, WorldPop data (https://www.worldpop.org/ (accessed on 6 June 2022)). In consideration of the inconformity between WorldPop data and those reported in the statistical yearbook, the WorldPop data employed a ratio coefficient method to perform the correction based on the population data of the statistical yearbook. The population weight of each grid can be calculated by the total population of the grid divided by the total population of the prefecture-level region covering the spatial grid. Then, the $PWGL{O}_3$ concentrations of each prefecture-level region may be calculated with the $GL{O}_3$ concentrations multiplying by the population weight and zonal accumulation.

2.1.2. Influencing Factors

The influencing factors of the $PWGL{O}_3$ concentrations mainly came from three aspects [39] illustrated by Figure 1: social economics, traffic development, and natural environment. Each category driver contains a certain number of proxy variables. The social economics involved 7 proxy variables: population density (PD), gross domestic product per capita (GDPPC), tertiary industry proportion (TIP), real estate investment per capita (REIPC), industrial electricity consumption per capita (IECPC), electricity consumption of residents per capita (ECRPC), and total social electricity consumption per capita (TSECPC). It should be noted that, in China, the tertiary industry mainly refers to the service industry, while the secondary industry refers to the mining, manufacturing, electricity (heating power, fuel gas, and water) production and supply, and construction industries. The traffic development included 3 proxy variables: bus ridership per capita (BRPC), number of taxis per capita (NTPC), and number of buses per capita (NBPC). The data for the above 10 proxy variables were collected from the China city statistical yearbook in corresponding years. The natural environment contained 9 proxy variables: urban green coverage rate (UGCR), urban construction area per capita (UCAPC), ratio of urban construction area (RUCA), green area per capita (GAPC), difference of summer average temperature from that in 2004 (DSAT), annual $P{M}_{2.5}$ concentrations (APMC), annual cumulative duration of sunshine (ACDS), annual average wind speed (AAWS), and annual average relative humidity (AARH). It should be noted that the DSAT was selected in order to underline the rising temperature. As previously mentioned, VOC and $N{O}_x$, the emissions of human activities that mainly include industry and traffic, are the common precursors of ozone. The three traffic factors in this paper include public traffic, bus, and taxi. Additionally, the factor pertaining to industry emissions lacked data availability. Consequently, the APMC, as a proxy variable of the emissions of the industry and private cars, was absorbed. Moreover, UGCR, UCAPC, RUCA, and GAPC, were collected from the China city statistical yearbook in corresponding years (https://data.cnki.net/trade/yearbook/Single/N2022040095?zcode=Z024 (accessed on 6 June 2022)). The raw grid data of annual average $P{M}_{2.5}$ and temperature with a spatial resolution of $10\mathrm{k}\mathrm{m}\times 10\mathrm{k}\mathrm{m}$ were downloaded from the National Earth System Science Data Center (http://www.geodata.cn/data/ (accessed on 6 June 2022)). The data of the three climate factors, ACDS, AAWS, and AARH, were collected from the China Meteorological Data Service Centre (http://data.cma.cn/en (accessed on 6 June 2022)).

2.2. Bayesian Multi-Stage Spatiotemporal Evolution Hierarchy Model (BMSSTEHM)

The BMSSTEHM presented by Li et al. [40] was employed to explore the spatiotemporal trends of the $PWGL{O}_3$ concentrations from 2005 to 2019 in mainland China. The BMSSTEHM integrated the BSTHM [41] and piecewise regression model [42] and can investigate non-linear local trends by capturing the self-adaptive turning points of local trends and considering spatial correlations and heterogeneity. The mathematical forms of the BMSSTEHM are as follows:

(2)$c_{it}~Normal({\mu }_{it},{\sigma }_c^2)$

(3)${\mu }_{it}=\alpha +S_i+\left(b_{0}t+v_t\right)+b_{1i}t+b_{2i}\left(t{-a}_{1i}\right)*G_{t,a_{1i}}+{\epsilon }_{i,t}$

(4)$G_{t,a_{1i}}=\frac{1}{1+\mathrm{e}\mathrm{x}\mathrm{p}(-\lambda (t-a_{1i}\left)\right)}$

(5)$a_{1\mathrm{i}}~\mathit{Uniform}(3,\mathrm{T}-4)$

(6)$k_{1i}={b_0+b}_{1i}$

(7)$k_{2i}={b_0+b}_{1i}+b_{2i}$

where $c_{it}$ represents the observed annual average $PWGL{O}_3$ concentrations of the i-th prefecture-level region at year $t$. ${\mu }_{it}$ and ${\sigma }_c^2$ represent the mean and variance of the likelihood distribution. $\mathsf{\alpha }$ is an intercept term. $S_i$ is the parameter of spatial fixed effect; $\left(b_{0}t+v_t\right)$ describes overall trends containing a linear $b_0$ and non-linear tendency $v_t$, whose prior distribution adapted Gaussian distribution; $b_{1i}$,$b_{2i}$ are local piecewise linear regression coefficients; $a_{1i}$ is the first turning point of the i-th prefecture-level region. If there is one turning point, the spatiotemporal can be divided into the first (early) and second (later) stages. The corresponding linear variation parameters are $k_{1i}$ and $k_{2i}$. $G_{t,a_{1i}}$ is the logistic function, as described in Equation (3). $\mathsf{\lambda }$ is a shape parameter that is generally assigned to a value greater than 10. $T$ is the number of time nodes ($T=15$ in this study). Linear regression is necessary only when there are more than two sample values, namely at least three sample values. $a_{1\mathrm{i}}$ is therefore limited to between 3 and $T-4$. Then, it can guarantee that there are at least three sample values on both linear sides. $a_{1\mathrm{i}}$ is assigned to the prior uniform distribution and can be self-adaptively estimated according to the sample data. The term ${\epsilon }_{it}$ is a Gaussian noise error whose prior distribution is assigned as a normal distribution $N(0,{\sigma }_{\epsilon }^2)$. The spatial relative magnitude of the $PWGL{O}_3$ concentrations of the i-th prefecture-level region in the two stages, denoted as ${SR}_{1i}$ and ${SR}_{2i}$, quantifying the $PWGL{O}_3$ polluted level relative to the global overall level can be estimated. ${SR}_{1i}$ and ${SR}_{2i}$ can be expressed as follows:

(8)${SR}_{1i}=\frac{{\widehat{\mu }}_{i,t_1^{*}}}{\frac{1}{M}{\sum }_{i=1}^M{\widehat{\mu }}_{i,t_1^{*}}}$

(9)${SR}_{2i}=\frac{{\widehat{\mu }}_{i,t_2^{*}}}{\frac{1}{M}{\sum }_{i=1}^M{\widehat{\mu }}_{i,t_2^{*}}}$

(10)${\widehat{\mu }}_{i,t_1^{*}}=\alpha +S_i+b_0{t}_1^{*}+b_{1i}{t}_1^{*}+{\epsilon }_1$

(11)${\widehat{\mu }}_{i,t_2^{*}}=\alpha +S_i+b_0{t}_2^{*}+b_{1i}{t}_2^{*}+b_{2i}\left(t_2^{*}{-a}_{1i}\right)*G_{t,a_{1i}}+{\epsilon }_2$

(12)$t_1^{*}=\frac{1+a_{1i}}{2}$

(13)$t_2^{*}=\frac{a_{1i}+T}{2}$

where $t_1^{\mathrm{*}}$ and $t_2^{*}$ represent the temporal centre points of the early and later stages, respectively. In this study, $M$ is the number of prefecture-level regions of mainland China, which is 338. The parameters, $S_i,$ $b_{1i}$,$b_{2i}$, and $a_{1\mathrm{i}}$, are all considered simultaneously spatial structured and unstructured effects by assigning the prior of the BYM model [43]. Therefore, the ‘Queen’ spatial adjacency structure sharing a common edge or a common point was adopted in this paper. ${\epsilon }_1$ and ${\epsilon }_2$ represent Gaussian error terms.

The Bayesian statistics calculations in this paper were implemented by WinBUGS 14.0 [44] based on the Markov chain Monte Carlo (MCMC) algorithm. The number of iterations was set to 200,000, 180,000 were the burn-in period, and 20,000 were the sampling number of the posterior distribution of parameters. The convergence of the Bayesian inferences in this study was monitored by standard autocorrelation plots and trace plots.

2.3. Bayesian Spatiotemporal LASSO Regression Model (BST-LASSO-RM)

For spatiotemporal process, the spatial endogenous and exogenous interaction effects should be considered [45,46,47]. Considering possible multicollinearity problems due to the no little number of explanatory variables, this study presented a BST-LASSO-RM based on the Bayesian LASSO regression model [48] to overcome the problem of multicollinearity, considering the spatial endogenous and exogenous interaction effects. The regression parameters of the BST-LASSO-RM were assigned to Laplace priors [49]. The BST-LASSO-RM can acquire a more stable estimation and automatically provide interval estimates for all parameters, including the error variance [48]. The mathematical structure of the BST-LASSO-RM can be expressed as follows:

(14)$y_{it}~Normal({\mu }_{it},{\sigma }_y^2)$

(15)${\mu }_{it}=\varphi +\rho {(W}_i{\mu }_{it})++{\sum }_{j=1}^k{\beta }_j{x}_{j,i,t}+{\sum }_{j=1}^k{\vartheta }_j{{(W}_{i}x}_{j,i,t})+{\xi }_{it}$

(16)$\widehat{\mathit{\beta }}=argmi{n}_{\beta }{\left(\mathit{\mu }-\mathit{\rho }\mathit{W}\mathit{\mu }-\mathit{\beta }\mathit{X}-\mathit{\vartheta }\mathit{W}\mathit{X}\right)}^T\left(\mathit{\mu }-\mathit{\rho }\mathit{W}\mathit{\mu }-\mathit{\beta }\mathit{X}-\mathit{\vartheta }\mathit{W}\mathit{X}\right)+{\lambda }_{\beta }{‖\mathit{\beta }‖}_1$

(17)$\widehat{\mathit{\rho }}=argmi{n}_{\rho }{\left(\mathit{\mu }-\mathit{\rho }\mathit{W}\mathit{\mu }-\mathit{\beta }\mathit{X}-\mathit{\vartheta }\mathit{W}\mathit{X}\right)}^T\left(\mathit{\mu }-\mathit{\rho }\mathit{W}\mathit{\mu }-\mathit{\beta }\mathit{X}-\mathit{\vartheta }\mathit{W}\mathit{X}\right)+{\lambda }_{\rho }{‖\mathit{\rho }‖}_1$

(18)$\widehat{\mathit{\vartheta }}=argmi{n}_{\vartheta }{\left(\mathit{\mu }-\mathit{\rho }\mathit{W}\mathit{\mu }-\mathit{\beta }\mathit{X}-\mathit{\vartheta }\mathit{W}\mathit{X}\right)}^T\left(\mathit{\mu }-\mathit{\rho }\mathit{W}\mathit{\mu }-\mathit{\beta }\mathit{X}-\mathit{\vartheta }\mathit{W}\mathit{X}\right)+{\lambda }_{\vartheta }{‖\mathit{\vartheta }‖}_1$

(19)$\mathit{\beta }|{\lambda }_{\beta },{\sigma }_y^2~{\prod }_{j=1}^k\frac{{\lambda }_{\beta }}{2{\sigma }_y}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{{\lambda }_{\beta }\left|{\beta }_j\right|}{{\sigma }_y})$

(20)$\mathit{\rho }|{\lambda }_{\rho },{\sigma }_y^2~\frac{{\lambda }_{\rho }}{2{\sigma }_y}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{{\lambda }_{\rho }\left|\rho \right|}{{\sigma }_y})$

(21)$\mathit{\vartheta }|{\lambda }_{\vartheta },{\sigma }_y^2~{\prod }_{j=1}^k\frac{{\lambda }_{\mathit{\vartheta }}}{2{\sigma }_y}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{{\lambda }_{\mathit{\vartheta }}\left|{\vartheta }_j\right|}{{\sigma }_y})$

where ${\mu }_{it}$ and ${\sigma }_y^2$ represent the mean and variance of the normal likelihood distribution of the observed $PWGL{O}_3$ concentrations. $\varphi$ is an intercept term. $W_i$ presents the spatial adjacent matrix of the $i$-th prefecture-level city and adopts Queen contiguity that includes common vertices and edges. $W_i{\mu }_{it}$ and ${W_{i}x}_{j,i,t}$ represent the spatial adjacent average $PWGL{O}_3$ concentrations and the j-th explaining variable around the i-th prefecture-level region in the t-th year. $\rho$ is the parameter of the spatial endogenous effect of the annual $PWGL{O}_3$ concentrations, $W_i{\mu }_{it}$. ${\beta }_j$ is the regression parameters of local influencing factors, $x_{j,i,t.}$, and ${\vartheta }_j$ is the parameter of spatial exogenous effects of spatial adjacent influencing factors, $W_i{x}_{j,i,t}$; $k$ is the number of the explanatory variables, ${\xi }_{it}$ represents the random error, $\widehat{\mathit{\beta }}$, $\widehat{\mathit{\rho }},$ and $\widehat{\mathit{\vartheta }}$ represent the estimate of the regression parameters $\mathit{\beta },$ $\mathit{\rho }$ and $\mathit{\vartheta }$. $\mathit{\mu }$ and $\mathit{X}$ are the matrices of ${\mu }_{it}$ and $x_{j,i,t}$, respectively. ${\lambda }_{\beta },$ ${\lambda }_{\rho },$ and ${\lambda }_{\vartheta }$ are the shrinkage parameters. The Bayesian calculations of the BST-LASSO-RM inferences were implemented mainly by the Python package, PyMC3.

2.4. A Spatiotemporal Propensity Score Matching (STPSM) Method

Similarly, to take into account the spatial endogenous and exogenous interaction effects, a STPSM method was presented based on the PSM method [50] in this study. This method can be adopted to recognise the causality between the $PWGL{O}_3$ concentrations, $\mathit{Y}$, and the spatial adjacent $PWGL{O}_3$ concentrations, $\mathit{W}\mathit{Y}$, the local and spatial adjacent influencing factors, $\mathit{X}$ and $\mathit{W}\mathit{X}$. Matched sets of treated and untreated groups that share a similar value of the propensity score can be formed by the STPSM, and then the average treatment effect for the treated, $E[Y\left(1\right)-Y(0\left)\right|Z=1]$ is allowed to be estimated. The propensity score, $Pr(Z_i=1|(W_i{y}_i,x_{i,1},\dots ,x_{i,k},{W_{i}x}_{i,1},\dots ,{W_{i}x}_{i,k}))$, is the probability of treatment assignment, $Z_i=1$, conditional on other observed covariates, ${\mathit{X}}_{\mathit{i}}$. The propensity score of one factor can be obtained by logistic regression as follows:

(22)$\mathit{Pr}\left(Z_i=1|(W_i{y}_i,x_{i,1},\dots ,x_{i,k},{W_{i}x}_{i,1},\dots ,{W_{i}x}_{i,k})\right)=\frac{1}{1+\exp \left(\rho W_i{y}_i+{\sum }_{k=1}^K(-{\mathit{\theta }}^T{\mathit{x}}_{i,k}-{\vartheta }^T{W}_i{x}_{i,k})\right)}$

Considering that all the variations in this study are continuous numeric variables, the sample will be divided into treated and untreated groups according to whether the target variate is greater than its median. STPSM is then employed to implement 1:1 matching, in which pairs of treated and untreated groups are formed with similar values of the propensity score.

2.5. A Spatiotemporal Causal Inference Method

This study proposed a method for recognizing spatiotemporal causality and cause–effect, integrating STPSM and BST-LASSO-RM. As shown in Figure 2, this method consists of three steps. First, the BST-LASSO-RM was employed to recognise the significant correlated, ${\mathsf{\Lambda }}^{\left(1\right)}$, or not significant correlated variables, ${\mathsf{\Lambda }}^{\left(0\right)}$, with the outcome variable, $\mathit{Y}$, from all independent variables, $\mathit{X}$, $\mathit{W}\mathit{Y}$, and $\mathit{W}\mathit{X}$.

(23)$y_{it}=\frac{Y_{it}-Y_{mean}}{Y_{std}}$

(24)$y_{i}t~Normal({\mu }_{it},{\sigma }_Y^2)$

(25)${\mu }_{it}=\rho W_i{\mu }_{it}+{\sum }_{j=1}^K{\beta }_j{x}_{j,i,t}+{\sum }_{j=1}^K{\vartheta }_j{W_{i}x}_{j,i,t}$

where $Y_i$ and $y_i$ are the observed value and $Z$ standardized score, $Y_{mean}$ and $Y_{std}$ are the corresponding mean and standard deviation. Specifically, the representations of the $W_i{\mu }_{it}$, ${\mathit{x}}_{j,i,t}$ and $W_i{x}_{j,i,t}$ are the same as the above. And the three type of variables can be identified as the significant correlated variable if the posterior probability of the corresponding regression coefficients, $P(\rho >0|Data)$ or $P(\rho <0|Data)$, ${P(\beta }_j>0\left|Data\right)$ or ${P(\beta }_j\left|Data\right)<0$, ${P(\vartheta }_j>0\left|Data\right)$ or ${P(\vartheta }_j\left|Data\right)<0$, are greater than $1-\alpha$, e.g., 0.95, and vice versa. Second, the STPSM method was used to recognise the causality of the significant correlated independent variables, ${\mathsf{\Lambda }}^{\left(1\right)}$ including $\mathit{W}\mathit{Y}$, $\mathit{X},$ and $\mathit{W}\mathit{X}$, to outcome variable, $\mathit{Y}$. Let us suppose there are $K(K\le k)$ significant correlated variables, ${\mathsf{\Lambda }}_j^{\left(1\right)}(j=\mathrm{1,2},\dots ,K)$, and then follow the steps below to recognise the causality of each of the ${\mathsf{\Lambda }}_j^{\left(1\right)}(j=\mathrm{1,2},\dots ,K)$:

• (i). Let ${\mathsf{\Lambda }}_{1,i,t}^{\left(1\right)}$ = 1 if ${\mathsf{\Lambda }}_{1,i,t}^{\left(1\right)}\ge Mean\left({\mathsf{\Lambda }}_{1,i,t}^{\left(1\right)}\right)$, ${\mathsf{\Lambda }}_1^{\left(1\right)}$ = 0 if ${\mathsf{\Lambda }}_{1,i,t}^{\left(1\right)}<Mean\left({\mathsf{\Lambda }}_{1,i,t}^{\left(1\right)}\right)$;

• (ii). Then, the observed outcome, $Y_{it}$, can be divided into two groups, they can be called treatment and control group, according to the value of ${\mathsf{\Lambda }}_{1,i,t}^{\left(1\right)}$, i.e., ${\mathsf{\Lambda }}_{1,i,t}^{\left(1\right)}=1$ and ${\mathsf{\Lambda }}_{1,i,t}^{\left(1\right)}=0$;

• (iii). Calculate propensity scores $P{S}_{{\mathsf{\Lambda }}_{1,i,t}^{\left(1\right)}=1}=f\left({\mathsf{\Lambda }}_{2,i,t}^{\left(1\right)},{\mathsf{\Lambda }}_{3,i,t}^{\left(1\right)},\dots ,{\mathsf{\Lambda }}_{k-1,i,t}^{\left(1\right)},{\mathsf{\Lambda }}_{k,i,t}^{\left(1\right)}\right)$;

• (iv). Using STPSM to match the sample data between the treatment and control groups;

• (v). Adopt t-test to identify the significant difference between the treatment and control groups and consequently to recognise the causality of ${\mathsf{\Lambda }}_1^{\left(1\right)}$ to $\mathit{Y}$;

• (vi). Repeat the above five steps, (i)–(v), to recognise the causality of other significant correlated independent variables, ${\mathsf{\Lambda }}_2^{\left(1\right)},\dots ,{\mathsf{\Lambda }}_M^{\left(1\right)}$.

The identified causal and not causal variables were denoted as ${\mathsf{\Lambda }}_j^{\left(11\right)}(j=\mathrm{1,2},\dots ,M\le K)$ and ${\mathsf{\Lambda }}_l^{\left(10\right)}(l=\mathrm{1,2},\dots ,K-M)$, it identifies that all significant correlated influencing factors are causal factors to explained variable, $Y$, if $K-M=0$. Through the above steps, possible three types of causal effects, spatial endogenous causal effect from $\mathit{W}\mathit{Y}$, spatial exogenous effect from $\mathit{W}\mathit{X}$, and local variational causal effect from $\mathit{X}$, can be identified. Third, again employing the BST-LASSO-RM to estimate the quantitative cause-effects, ${\mathit{\delta }}_{{\mathsf{\Lambda }}_j^{\left(11\right)}}$:

(26)$\mathit{Y}=\psi +{\sum }_j^K{\mathit{\delta }}_{{\mathsf{\Lambda }}_j^{\left(11\right)}}{\mathsf{\Lambda }}_j^{\left(11\right)}+\mathit{\epsilon }$

3. Results

3.1. Descriptive Statistics

Overall, Chinese ground-level $O_3$ concentrations maintain a slight downward trend from 2005 to 2014, but an increasing tendency arose after 2014. The median of ground-level $O_3$ annual concentrations decreased from 82.73 in 2005 to 81.57 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$ in 2014, but increased to 87.84 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$ in 2019. The $PWGL{O}_3$ annual concentrations maintained relatively steady from 2005 (81.80 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$) to 2014 (81.06 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$), an increasing trend occurred from 2014 to 2019 (89.92 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$). The heterogeneity of the $PWGL{O}_3$ concentrations at the prefectural level in China showed an overall rising trend. The corresponding coefficient of variation (standard deviation divided by mean) increased from 14.16% in 2005 to 15.60% in 2014, and to 17.48% in 2019. Figure S1 shows the temporal series of the spatial distribution of $PWGL{O}_3$ annual concentrations in mainland China from 2005 to 2019. This reveals that the spatial pattern of the $PWGL{O}_3$ annual concentrations in mainland China has changed since 2014. The number of prefecture-level regions, with the annual $PWGL{O}_3$ concentrations greater than 80 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$, increased from 172 in 2014 to 244 in 2019. The $PWGL{O}_3$ concentrations in the North China Plain region have become even more serious since 2014. In other words, different spatial structures of the $PWGL{O}_3$ concentrations were formed during 2005–2014 and 2015–2019.

3.2. The Spatial Pattern and Its Transformation of the $PWGL{O}_3$

The spatiotemporal process can be generally divided into two phases called early and late stages. The overall spatial distributions of the $PWGL{O}_3$ concentrations in the two phases can be identified by the BMSSTEHM (Figure 3). The results show that the overall spatial pattern of the $PWGL{O}_3$ concentrations transformed from the early to late stage. Specifically, the spatial distribution of the $PWGL{O}_3$ concentrations in the early stage is disperser than that in the late stage. The high $PWGL{O}_3$ polluted regions transferred from the eastern and southern coastal areas in the early stage into the northern and North China Plain regions, including Hebei, Shanxi, Shandong, Henan, Jiangsu, and Anhui, in the late stage. The composition of the number of prefecture-level regions with the five classifications of the $PWGL{O}_3$ levels changed between the two stages. The numbers of the high $PWGL{O}_3$ polluted prefecture-level regions with spatial relative magnitude (SRM) of 1.05~1.20 and greater than 1.20 changed from 110 (32.5%) and 13 (3.8%) in the early stage to 94 (27.8%) and 68 (20.1%) in the late stage. The number of the medium $PWGL{O}_3$ polluted prefecture-level regions decreased from 96 (28.4%) in the early stage to 64 (18.9%) in the late stage. The number of the lightly $PWGL{O}_3$ polluted prefecture-level regions with SRM of less than 0.95 remained stable during the study period, 2005–2019.

Furthermore, the SRM of $PWGL{O}_3$ concentrations over the southeast coast areas of China decreased from 1.05~1.20 in the early stage to less than 1.05 in the late stage, while that over the North China Plain regions increased to greater than 1.20. The $PWGL{O}_3$ concentrations SRM over the northern regions, such as Shanxi, Inner Mongolia, Gansu and Ningxia, increased from less than 1.05 in the early stage to 1.05~1.20 and greater than 1.20 in the late stage.

3.3. The Local Trends of the $PWGL{O}_3$ Concentrations at the Sub-Provincial Level

The annual change in the $PWGL{O}_3$ concentrations in the two stages at the prefecture level can also be estimated by the parameters, $b_{1i}$ and $b_{1i}+b_{2i}$, of the BMSSTEHM. Research shows that the $GL{O}_3$ production is a photochemical oxidation process determined highly by meteorological conditions and precursors, such as VOC and $N{O}_x$. Consequently, variations in the annual mean $PWGL{O}_3$ concentrations could be modulated by many factors, especially climate factors, such as air temperature and sun radiation, which can elevate or lower annual mean $GL{O}_3$ concentrations. To validate the reliability of the annual temporal trends estimated by the BMSSTEHM, the scatters and BMSSTEHM fitted polylines of the annual $PWGL{O}_3$ concentrations in the example of 12 prefecture-level regions illustrated in Figure 4. The 12 prefecture-level regions were randomly selected while ensuring coverage in all directions of China. The results show that the fitted polylines can capture different linear trends in the two stages, although the annual $PWGL{O}_3$ concentrations in some years were saltatory or abnormal, possibly due to unstable climate factors. It should also be noted that the parameters of the BMSSTEHM, $b_{1i}$ and $b_{1i}+b_{2i}$, measured the average annual change within the corresponding two stages; thus, the mutations in some years had little influence on the two parameters.

Figure 5 illustrates the local trends of the $PWGL{O}_3$ concentrations in the early stages. Generally, most prefecture-level regions (72.2%) experienced a decreased local trend of the $PWGL{O}_3$ concentrations in the early stages. The corresponding median, lower, and upper quartiles of the annual change are −0.36, −0.53, and −0.19 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$ per year. The other 27.8% of the prefecture-level regions are distributed dispersedly across China’s mainland. The median, lower, and upper quartiles of the annual increasing changes are 0.20, 0.08, and 0.64 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$ per year. Some sub-provincial regions located in Guangxi, Guangdong, and Tibet possessed high $PWGL{O}_3$ polluted levels and increased local trends concurrently in the early stage.

In the late stage, the spatial pattern of the local trends of the $PWGL{O}_3$ concentrations formed a distinct spatial structure (Figure 6). The $PWGL{O}_3$ concentrations in the most sub-provincial areas (79.3%) showed an increased trend. In particular, some cities (dark red coloured areas in Figure 6) of central China, such as Ningxia, Shanxi, Henan, Anhui, and south Gansu, experienced a strong increased trend of greater than 4.00 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$ per year. The median, lower, and upper quartiles of all the prefecture-level regions with increasing local trends are 2.82, 1.19, and 4.38 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$ per year. The local trends of the areas located in Tibet, Xinjiang, Guangxi, and Guangdong converted from upward in the early stage to downward trend in the late stage.

3.4. Drivers of the $PWGL{O}_3$ Concentrations

First, the BST-LASSO-RM presented in this paper was employed to assess the significance of the correlation between the $PWGL{O}_3$ concentrations and the 19 local independent variables (PD, GDPPC, TIP, REIPC, IECPC, ECRPC, TSECPC, BRPC, NTPC, NBPC, UGCR, GAPC, UCAPC, RUCA, DSAT, APMC, ACDS, AAWS, and AARH) and 20 corresponding spatial adjacent variables (W·$\mathrm{P}\mathrm{W}\mathrm{G}\mathrm{L}{\mathrm{O}}_3$ and the other 19 spatial adjacent variables, e.g., W·PD, W·GDPPC, etc.) in the two stages, 2005–2014 and 2015–2019. According to the results estimated by the BST-LASSO-RM, 13 local variables (GDPPC, TIP, ECRPC, TSECPC, BRPC, NTPC, NBPC, GAPC, DSAT, APMC, ACDS, AAWS, and AARH) and 4 spatial adjacent variables (W·$\mathrm{P}\mathrm{W}\mathrm{G}\mathrm{L}{\mathrm{O}}_3,$ W·GDPPC, W·TSECPC, and W·AAWS) were significantly associated with the $PWGL{O}_3$ concentrations in 2005–2014. In the period of 2015–2019, the significant influence pattern changed to 13 local variables (GDPPC, TIP, ECRPC, TSECPC, BRPC, NTPC, NBPC, UGCR, DSAT, APMC, ACDS, AAWS, and AARH) and 7 spatial adjacent variables (W·$\mathrm{P}\mathrm{W}\mathrm{G}\mathrm{L}{\mathrm{O}}_3,$ W·GDPPC, W·TIP, W·ECRPC, W·TSECPC, W·APMC, W·AAWS), which are significant.

Secondly, recognising the causality of the significant influence factors on the $PWGL{O}_3$ concentrations. As mentioned in Section 2, the STPSM method was adopted to recognise causality. Considering that all the variables are numeric types for each significant influence variable, the prefecture-level regions can be bounded by the median to be divided into the treatment group and the control group. Then, the prefecture-level regions in the treatment group and control group were matched by the STPSM. Subsequently, the relative risk (RR) and risk difference (RD) of each selected significant influence variable on the explanatory variable, the $PWGL{O}_3$ concentrations, can be estimated (Table 1). Furthermore, the RR and RD can be calculated by the ratio and difference of the average $PWGL{O}_3$ concentrations in the matched treatment group and control group. The results show that most of the selected significant influence variables possessed statistically significant (p < 0.05) causality with the $PWGL{O}_3$ concentrations in 2005–2014 and 2015–2019. Precisely, W·$\mathrm{P}\mathrm{W}\mathrm{G}\mathrm{L}{\mathrm{O}}_3,$ GDPPC, W·GDPPC, TSECPC, W·TSECPC, GAPC, DSAT, APMC, ACDS, and W·AAWS are risk factors, and BRPC, NTPC, NBPC, and ARRH are protective factors in 2005–2014. To the late stage, 2015–2019, the risk factors increased from 10 to 16 for W·$\mathrm{P}\mathrm{W}\mathrm{G}\mathrm{L}{\mathrm{O}}_3,$ GDPPC, W·GDPPC, TIP, W·TIP, ECRPC, W·ECRPC, TSECPC, W·TSECPC, NTPC, UGCR, DSAT, APMC, W·APMC, ACDS, AAWS, W·AAWS, and AARH, and the protective factors decreased from 4 to 3, BRPC, AAWS, and AARH.

Figure 7 shows the normalised regression coefficients of the BST-LASSO-RM between the $PWGL{O}_3$ concentrations and the driving factors in the two stages. The results indicate that the structure of the drivers of the $PWGL{O}_3$ concentrations has changed from the early stage to the late stage. Concretely, the eight factors, GDPPC, TSECPC, BRPC, NTPC, DSAT, ACDS, ARRH, and three spatial adjacent factors, W·$\mathrm{P}\mathrm{W}\mathrm{G}\mathrm{L}{\mathrm{O}}_3$, W·GDPPC, W·TSECPC, had driving effects in both stages. However, the driving direction of the NTPC factor changed from negative to positive. This is possibly due to the fact that the use efficiency of taxis was reduced from the early to the late stage due to the increase in private cars. Moreover, GAPC and NBPC had driving effects in 2005–2014, excluding 2015–2019; however, TIP, W·TIP, ECRPC, W·ECRPC, UGCR, W·APMC, AAWS, and W·AAWS were on the contrary. A finding is that GAPC has a positive driving effect in only one stage. The reason for this could be that GAPC decreased from the early stage to the late stage. The negative driving effects in 2005–2014 for BRPC, NTPC, and NBPC, and 2015–2019 for BRPC, indicate that the development of public transportation, especially bus ridership, can alleviate $PWGL{O}_3$ concentrations. Six local factors, BRPC, APMC, TSECPC, GDPPC, NBPC, and GAPC, were the major drivers in 2005–2014, and the corresponding explanatory powers, the percentages of the absolute value of the normalised regression coefficients, were 13.56%, 11.95%, 9.82%, 8.72%, 8.46%, and 8.27%. Five local factors and one spatial adjacent factor, APMC, ACDS, TSECPC, BRPC, AARH, and W·APMC, were the major drivers with explanatory powers of 17.57%, 9.75%, 9.71%, 7.05%, 6.28%, and 8.94%, in 2015–2019.

To estimate the causal effects of the driving factors on the $PWGL{O}_3$ concentrations across China in the two stages, the non-normalised regression results of the BST-LASSO-RM were estimated (Table 2). The results indicate that the driving effects of the three spatial adjacent drivers, W·$\mathrm{P}\mathrm{W}\mathrm{G}\mathrm{L}{\mathrm{O}}_3$, W·GDPPC, and W·TSECPC, increased from 0.76 (0.18, 1.17), 0.12 (0.05, 0.51), and 0.62 (0.22, 1.08), in 2005–2014 to 1.86 (2.55, 2.19), 0.33 (0.12, 0.56), and 1.94 (1.01, 2.30), in 2015–2019. The driving effects of GDPPC and BRPC decreased from 0.45 (0.23, 0.76) and −3.11 (−5.15, −0.99), in 2005–2014 to 0.19 (0.26, 0.39) and −0.96 (−1.31, −0.58), in 2015–2019. The driving effects of TSECPC, DSAT, APMC, ACDS, and ARRH, increased from 0.91 (0.42, 1.54), 0.08 (0.03, 0.12), 0.03 (0.01, 0.06), and −0.07 (−0.12, −0.03) in 2005–2014 to 2.34 (1.67, 3.18), 0.16 (0.09, 0.21), 0.12 (0.06, 0.22), and −0.12 (−0.36, −0.10) in 2015–2019. The driving effect of NTPC decreased, but the driving direction reversed. The UGCR and AAWS have positive driving effects, 0.04 (0.01, 0.06) and 0.08 (0.04, 0.18), on the $PWGL{O}_3$ concentrations in 2015–2019, excluding 2005–2014.

3.5. Implications of the Results

The results of the transitions of the overall spatial pattern and local trends of the annual $PWGL{O}_3$ concentrations can clarify the spatiotemporal evolution mechanism of the $PWGL{O}_3$ concentrations in China. This indicates that the emphatically concerned regions should focus on the areas with higher levels and increasing trends of $PWGL{O}_3$ concentrations (i.e., the red and dark red areas in Figure 3B and Figure 5). The identification of the potential drivers regarding the prefecture level of the $PWGL{O}_3$ concentrations in China can help policy-makers design more effective policies to control $PWGL{O}_3$ concentrations. The annual $PWGL{O}_3$ concentrations can be reduced by decreasing the factors with positive driving effects and increasing the factors with negative driving effects. Specifically, although local and spatial adjacent social economic factors (GDPPC, W·GDPPC, TIP, W·TIP, ECRPC, W·ECRPC, TSECPC, and W·TSECPC) have positive effects, these factors cannot be easily decreased due to the development of the economy. Nevertheless, the annual $PWGL{O}_3$ concentrations can be lowered by decreasing NTPC and increasing BRPC simultaneously. The increase in UGCR can result in the growth of the annual $PWGL{O}_3$ concentrations, and green coverage may provide many human health benefits. This finding suggests that the vegetation, with more $VOC$ emissions as precursors of ozone, should be cut down in the urban area. The results estimated by the spatiotemporal causal inference method show that the spatial adjacent $PWGL{O}_3$, i.e., W·$PWGL{O}_3$, has positive driving effects in the two stages, 2005–2014 and 2015–2019, this states that the $PWGL{O}_3$ concentrations have a significant spatial causal spillover effect. Namely, the local $PWGL{O}_3$ concentrations in a prefecture-level region may also be aggravated by its spatial adjacent $PWGL{O}_3$ concentrations. Meanwhile, considering that the AMPC and W·APMPC have positive driving effects, this indicates that surface ozone pollution may be synergistically controlled with local and circumjacent $P{M}_{2.5}$ pollution, i.e., synergistic governance. The climate drivers, such as temperature, sunshine, wind speed, and relative humidity, can also not be changed; however, the related results may help us understand the variation of $PWGL{O}_3$ concentrations in the background of global climate change.

4. Discussion

This paper focused on the problem of population exposure to surface ozone concentrations, measured by the annual $PWGL{O}_3$ concentrations, at a prefecture-level scale across China mainland during 2005–2019. In this study, an advanced method called BMSSTEHM was used to carefully detect the space–time trends of the annual $PWGL{O}_3$ concentrations. Furthermore, to explore the drivers of the $PWGL{O}_3$ concentrations, this study transcended the constraints of traditional correlation analyses by presenting a spatiotemporal causal inference method that discerns causality relationship and estimates the cause—effect. This method combines STPSM and BST-LASSO-RM, both of which are also proposed in this paper. The methods in this paper provide a methodological reference for other similar studies. And the spatiotemporal causal inference method can be applied potentially in analysing remote sensing data. Furthermore, the results of the paper can help in the creation of diverse regional policies to reduce population exposure to $GL{O}_3$ concentrations.

The spatial distribution of $PWGL{O}_3$ concentrations in mainland China has shifted since 2015. It is important to carry out in-depth research on why the spatial pattern of $PWGL{O}_3$ concentrations on China’s mainland has changed. The spatial structure of the local trends of the annual $PWGL{O}_3$ concentrations at prefecture-level in China’s mainland showed differences in the two stages. The results indicated that the annual $PWGL{O}_3$ concentrations in the early stage maintained a steady state, i.e., most prefecture-level regions experienced a low local trend. But in the later period, the regions with high $PWGL{O}_3$ concentrations level experienced simultaneously high local trends. This means that the $PWGL{O}_3$ concentrations of the hot spot regions will become more severe. Additionally, the areas with high levels of $PWGL{O}_3$ concentrations have the tendency to expand outwards.

The transformation reasons for the static spatial pattern and local trends of the annual $PWGL{O}_3$ concentrations are very intricate. Generally, the reasons can be classified as anthropogenic and meteorological influences. Some researchers [51,52,53] have concluded that meteorological factors were not negligible but anthropogenic factors were dominant. This study argues that the economic development and remaining higher $P{M}_{2.5}$ concentrations jointly lead to the transformation of spatiotemporal trends of Chinese $PWGL{O}_3$ concentrations. According to Wei et al.’s newest study [54], the $P{M}_{2.5}$ concentrations across China decreased significantly from 2015–2020, but those in North China still remained at a higher level; in contrast, those in South China reached a lower level. In particular, the Pearl River Delta region of China has carried out pioneering pollution control measures since 2014 (https://www.mee.gov.cn/ywgz/xcjy/shxc/201909/t20190909_733028.shtml (accessed on 6 January 2023)), been “withdrawn” from the three key regions of air pollution prevention and control in China. The results of this study indicated that the $P{M}_{2.5}$ as a precursor of ozone has a maximal driving effect on $PWGL{O}_3$ in the late stage. Moreover, the rapid economic development of North China in the late stage may be also another reason.

The formation mechanism of the $PWGL{O}_3$ concentrations is very complex. Our study has preliminarily explored the drivers of the annual $PWGL{O}_3$ concentrations at the prefecture-level scale. The driving pattern of the annual $PWGL{O}_3$ concentrations in China’s mainland showed different features in the two phases. The four local and spatial adjacent economic factors, GDPPC, W·GDPPC, TIP, W·TIP, ECRPC, W·ECRPC, TSECPC, and W·TSECPC, all exhibited positive driving effects on the $PWGL{O}_3$ concentrations in the late stage. This indicates that the $P{M}_{2.5}$ concentrations were reduced by the shrinkage of the secondary industry due to environmental regulations; however, the development of the tertiary industry led to increasing $GL{O}_3$ concentrations. Moreover, not only the local economic factor can drive the $PWGL{O}_3$ concentrations, but also of the corresponding surrounding areas also can. This implies that the process of adjusting the policies aimed to alleviating surface ozone pollution must consider regional coordination.

In terms of traffic development factors, increasing bus ridership can promote the reduction in the $PWGL{O}_3$ concentrations in the two stages, and the driving effects of the BRPC in the early and later stages came close. It is interesting that the increase in NTPC and NBPC can lower the $PWGL{O}_3$ concentrations in the early stage; however, in the later stage, NBPC had no significant driving effect and the NTPC’s driving direction changed from negative to positive. The reason for the above phenomenon might be that travel demand within the city has increased rapidly due to the development of the tertiary industry. But the transport efficiency of the taxi is analogous to that of the private car, and less than of the bus. The low transport efficiency of the taxi may lead to high per capita exhaust gas emission, a precursor of ozone pollution. This suggests that public traffic should be energetically developed to improve public transport volumes and then reduce the $GL{O}_3$ concentrations.

This study found that GAPC and UGCR had positive driving effects in the early and later stages, respectively. This result is consistent with the conclusion of related environmental chemistry research [55]. Furthermore, the $GL{O}_3$ is mainly generated by photochemical reactions of $N{O}_x$ and $VOCs$ that were released from anthropogenic activities [56], while the $VOCs$ are an emission from vegetation [55]. Hence, green areas or green covers become risk factors for the $GL{O}_3$ concentrations. Our study quantified the relationship between the annual $PWGL{O}_3$ concentrations and GAPC and UGCR. Eucalyptus, Populus, Cunninghamia, etc., are plants with high VOC emissions; however, thickets and grasslands have low VOC emissions [57]. This result suggests that plants with high VOC emissions should be decreased in the urban green coverage.

For the factors of traffic development and green cover, the corresponding driving effects only occurred in the local area. That is, the factors of traffic development and green cover in the spatial adjacent areas have not driving effects on the local area’s $PWGL{O}_3$ concentrations. This finding provided important evidence or reference to control surface ozone pollution by changing the local factors of the traffic development and green cover.

From the view of ecological studies, this paper also demonstrated the conclusion that temperature and $P{M}_{2.5}$ are important drivers on the $GL{O}_3$ concentrations. We evaluated the quantitative driving effects of the local DSAT and APMC, and spatial adjacent $P{M}_{2.5}$, W·APMC, on annual $PWGL{O}_3$ concentrations at the prefectural scale, which is different from an environmental chemistry study. This can provide valuable evidence for policy-making for mitigating $PWGL{O}_3$ concentrations in China and other global regions.

This paper contributes by providing not only the concrete statistical results of the spatiotemporal trends of the $PWGL{O}_3$ concentrations at a sub-provincial scale, but also a frame for recognising causality and estimating causal effects for the situation with multi-factors, and investigating the drivers of the $PWGL{O}_3$ concentrations in mainland China. However, this study also has some limitations. First, the study period was not comprehensive enough, and the statistical results would be more robust if the study period was long enough. Second, the spatial scale was adopted with sub-province, but this can still be refined, such as the county scale. The above limitations will be emphatically focused on.

5. Conclusions

First, this study proposed a spatiotemporal causal inference method combining STPSM and BST-LASSO-RM, both of which are also proposed in this paper, and this method can be applied potentially in analyzing remote sensing spatiotemporal data. Second, the spatial distribution of $PWGL{O}_3$ concentrations in mainland China has shifted since 2015. Third, while early-stage annual $PWGL{O}_3$ concentrations remained stable across most prefecture-level regions in China’s mainland, the later stage the regions with high $PWGL{O}_3$ concentrations also exhibited high local trends. Fourth, a shift in the drivers of $PWGL{O}_3$ concentrations occurred from the early to late stages, the factors like NTPC changing their driving direction. Notably, the development of public transportation, especially bus ridership, can mitigate $PWGL{O}_3$ concentrations, with six local factors being the primary drivers in 2005–2014 and five local factors, along with one spatial adjacent factor, dominating in 2015–2019. Fifth, the annual $PWGL{O}_3$ concentrations can be lowered by decreasing NTPC and increasing BRPC simultaneously. Synergistic governance with local and adjacent $P{M}_{2.5}$ pollution can help control surface ozone pollution in the context of global climate change, while economic factors and climate drivers are challenging to modify.

Author Contributions

Conceptualization, J.L. and Z.R.; Methodology, J.L.; Software, J.L., Y.Y. (Yiming Yu), H.A., X.Y. and Y.Y. (Yixue Yang); Validation, J.X.; Formal analysis, J.L. and J.X.; Investigation, J.X.; Resources, J.W.; Data curation, Z.R., Y.Y. (Yiming Yu), H.A., X.Y. and Y.Y. (Yixue Yang); Writing—original draft, J.L.; Writing—review & editing, J.W. and Z.R.; Supervision, J.W. and Z.R.; Project administration, J.L.; Funding acquisition, J.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

Our deepest gratitude is expressed to our anonymous reviewers and editors for their careful work and constructive suggestions that have helped improve our paper substantially. We thank the support from the Save 2050 Programme jointly sponsored by Swarma Club and X-Order.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. The schematic diagram between outcome variable (dependent variable) and determinants and its proxy variables.

Figure 2. Flow diagram of the spatiotemporal causal inference framework.

Figure 3. The overall spatial pattern of the [Forumla omitted. See PDF.] concentrations in the early stage (A) and late-stage (B) and the transformation (C) of the composition of the five classes of [Forumla omitted. See PDF.] levels between the two stages.

Figure 4. The scatters and Bayesian multi-stage spatiotemporal evolution hierarchy model fitted polylines of the annual population-weighted ozone ([Forumla omitted. See PDF.]) concentrations in the example of 12 prefecture-level regions in China from 2005 to 2019.

Figure 5. The local annual change of the [Forumla omitted. See PDF.] concentrations in China at the sub-provincial scale in the early stage.

Figure 6. The local annual change in the [Forumla omitted. See PDF.] concentrations of China at a sub-provincial scale in the late stage.

Figure 7. Normalised regression results of the Bayesian LASSO regression model between the [Forumla omitted. See PDF.] concentrations and the significant driving factors in the two stages.

Table 1

The estimated relative risk (RR) and risk difference (RD) of the significant influence factors in the two periods, 2005–2014 and 2015–2019, based on the STPSM eliminating confounding effects.

Variables	2005–2014			2015–2019
	RR	RD (μg/m3)	p Value	RR	RD (μg/m3)	p Value
$\mathrm{W}·\mathrm{P}\mathrm{W}\mathrm{G}\mathrm{L}{\mathrm{O}}_3$	1.011	1.18	0.019	1.112	2.08	0.002
$\mathrm{G}\mathrm{D}\mathrm{P}\mathrm{P}\mathrm{C}$	1.025	2.03	0.017	1.023	2.15	0.018
$\mathrm{W}·\mathrm{G}\mathrm{D}\mathrm{P}\mathrm{P}\mathrm{C}$	1.016	1.19	0.011	1.015	2.08	0.008
TIP	1.001	0.52	0.289	1.044	3.68	<0.001
$\mathrm{W}·\mathrm{T}\mathrm{I}\mathrm{P}$	/	/	/	1.026	2.66	0.001
$\mathrm{E}\mathrm{C}\mathrm{R}\mathrm{P}\mathrm{C}$	1.025	1.87	0.158	1.042	3.81	0.048
$\mathrm{W}·\mathrm{E}\mathrm{C}\mathrm{R}\mathrm{P}\mathrm{C}$	/	/	/	1.021	2.78	<0.001
$\mathrm{T}\mathrm{S}\mathrm{E}\mathrm{C}\mathrm{P}\mathrm{C}$	1.061	5.07	<0.001	1.106	8.96	<0.001
$\mathrm{W}·\mathrm{T}\mathrm{S}\mathrm{E}\mathrm{C}\mathrm{P}\mathrm{C}$	1.022	3.57	0.006	1.058	6.61	0.002
$\mathrm{B}\mathrm{R}\mathrm{P}\mathrm{C}$	0.973	−3.20	0.001	0.957	−4.97	<0.001
$\mathrm{N}\mathrm{T}\mathrm{P}\mathrm{C}$	0.962	−5.19	<0.001	1.048	4.30	<0.001
$\mathrm{N}\mathrm{B}\mathrm{P}\mathrm{C}$	0.958	−3.20	<0.001	1.002	0.49	0.281
$\mathrm{U}\mathrm{G}\mathrm{C}\mathrm{R}$	/	/	/	1.044	3.68	<0.001
$\mathrm{G}\mathrm{A}\mathrm{P}\mathrm{C}$	1.036	2.75	0.009	/	/	/
$\mathrm{D}\mathrm{S}\mathrm{A}\mathrm{T}$	1.024	1.88	0.015	1.008	0.76	0.016
$\mathrm{A}\mathrm{P}\mathrm{M}\mathrm{C}$	1.043	3.34	0.001	1.150	11.75	<0.001
$\mathrm{W}·\mathrm{A}\mathrm{P}\mathrm{M}\mathrm{C}$	/	/	/	1.089	6.87	0.001
$\mathrm{A}\mathrm{C}\mathrm{D}\mathrm{S}$	1.006	1.88	0.011	1.076	2.88	<0.001
$\mathrm{A}\mathrm{A}\mathrm{W}\mathrm{S}$	/	/	/	0.951	−2.06	0.047
W$·$AAWS	1.002	1.061	0.125	1.008	1.94	0.033
AARH	0.991	−0.85	0.038	0.950	−2.25	0.002

Table 2

Non-normalised regression results of the Bayesian LASSO regression model between the $PWGL{O}_3$ concentrations and the significant driving factors in the two stages.

Drivers	2005–2014	2015–2019
	Mean (2.5%th, 97.5%th)	Mean (2.5%th, 97.5%th)
$\mathrm{W}\text{·}\mathrm{P}\mathrm{W}\mathrm{G}\mathrm{L}{\mathrm{O}}_3$$(\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$)	0.76 (0.18, 1.17)	1.86 (0.55, 2.19)
$\mathrm{GDPPC}({10}^3$ China Yuan)	0.45 (0.23, 0.76)	0.19 (0.26, 0.39)
$\mathrm{W}\text{·}\mathrm{GDPPC}({10}^3$ China Yuan)	0.12 (0.05, 0.51)	0.33 (0.12, 0.56)
$\mathrm{TIP}(\%$)	/	0.12 (0.07, 0.17)
W·TIP (%)	/	0.08 (0.03, 0.14)
$\mathrm{ECRPC}(\mathrm{k}\mathrm{W}\mathrm{h}$)	/	0.81 (0.07, 2.21)
$\mathrm{W}\text{·}\mathrm{ECRPC}(\mathrm{k}\mathrm{W}\mathrm{h}$)	/	0.62 (0.37, 0.99)
$\mathrm{TSECPC}({10}^2\mathrm{k}\mathrm{W}\mathrm{h}$)	0.91 (0.42, 1.54)	2.34 (1.67, 3.18)
$\mathrm{W}\text{·}\mathrm{TSECPC}({10}^2\mathrm{k}\mathrm{W}\mathrm{h}$)	0.62 (0.22, 1.08)	1.94 (1.01, 2.30)
$\mathrm{BRPC}({10}^2$ Passengers)	−3.11 (−5.15, −0.99)	−0.96 (−1.31, −0.58)
NTPC (Vehicles)	−0.13 (−0.30, −0.03)	0.08 (0.01, 0.16)
NBPC (Vehicles)	−0.14 (−0.30, −0.02)	/
UGCR (%)	/	0.04 (0.01, 0.06)
$\mathrm{GAPC}({\mathrm{m}}^2$)	0.61 (0.15, 1.11)	/
$\mathrm{DSAT}(\mathrm{℃}$)	0.08 (0.03, 0.12)	0.16 (0.09, 0.21)
$\mathrm{APMC}(\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$)	0.31 (0.20, 0.44)	0.65 (0.55, 0.74)
$\mathrm{W}\text{·}\mathrm{APMC}(\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$)	/	0.29 (0.09, 0.80)
$\mathrm{ACDS}(10\mathrm{D}\mathrm{a}\mathrm{y}\mathrm{s}$)	0.03 (0.01, 0.06)	0.12 (0.06, 0.22)
$\mathrm{AAWS}(\mathrm{m}/\mathrm{s}$)	/	0.08 (0.04, 0.18)
$\mathrm{W}\text{·}\mathrm{AAWS}(\mathrm{m}/\mathrm{s}$)	/	0.04 (0.02, 0.07)
$\mathrm{ARRH}(\%$)	−0.07 (−0.12, −0.03)	−0.21 (−0.36, −0.10)

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs15194871/s1, Figure S1: Temporal series of spatial distribution of ground-level ozone concentrations ($\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$) of mainland China from 2005 to 2019.

References

1. Agathokleous, E.; Feng, Z.; Oksanen, E.; Sicard, P.; Wang, Q.; Saitanis, C.J.; Araminiene, V.; Blande, J.D.; Hayes, F.; Calatayud, V. Ozone affects plant, insect, and soil microbial communities: A threat to terrestrial ecosystems and biodiversity. Sci. Adv.; 2020; 6, eabc1176. [DOI: https://dx.doi.org/10.1126/sciadv.abc1176] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32851188]

2. Lefohn, A.S.; Malley, C.S.; Smith, L.; Wells, B.; Hazucha, M.; Simon, H.; Naik, V.; Mills, G.; Schultz, M.G.; Paoletti, E. Tropospheric ozone assessment report: Global ozone metrics for climate change, human health, and crop/ecosystem research. Elem. Sci. Anthr.; 2018; 6, 27.

3. Sicard, P.; Anav, A.; De Marco, A.; Paoletti, E. Projected global ground-level ozone impacts on vegetation under different emission and climate scenarios. Atmos. Chem. Phys.; 2017; 17, pp. 12177-12196.

4. Stocker, T. Climate Change 2013: The Physical Science Basis: Working Group I Contribution to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2014.

5. Sicard, P. Ground-level ozone over time: An observation-based global overview. Curr. Opin. Environ. Sci. Health; 2021; 19, 100226.

6. Feng, L.; Liao, W. Legislation, plans, and policies for prevention and control of air pollution in China: Achievements, challenges, and improvements. J. Clean. Prod.; 2016; 112, pp. 1549-1558.

7. Xue, W.; Wang, L.; Yang, Z.; Xiong, Z.; Li, X.; Xu, Q.; Cai, Z. Can clean heating effectively alleviate air pollution: An empirical study based on the plan for cleaner winter heating in northern China. Appl. Energy; 2023; 351, 121923. [DOI: https://dx.doi.org/10.1016/j.apenergy.2023.121923]

8. Wang, Z.; Lv, J.; Tan, Y.; Guo, M.; Gu, Y.; Xu, S.; Zhou, Y. Temporospatial variations and Spearman correlation analysis of ozone concentrations to nitrogen dioxide, sulfur dioxide, particulate matters and carbon monoxide in ambient air, China. Atmos. Pollut. Res.; 2019; 10, pp. 1203-1210. [DOI: https://dx.doi.org/10.1016/j.apr.2019.02.003]

9. Cheng, L.; Wang, S.; Gong, Z.; Li, H.; Yang, Q.; Wang, Y. Regionalization based on spatial and seasonal variation in ground-level ozone concentrations across China. J. Environ. Sci.; 2018; 67, pp. 179-190. [DOI: https://dx.doi.org/10.1016/j.jes.2017.08.011]

10. Lyu, Y.; Wu, Z.; Wu, H.; Pang, X.; Qin, K.; Wang, B.; Ding, S.; Chen, D.; Chen, J. Tracking long-term population exposure risks to PM2. 5 and ozone in urban agglomerations of China 2015–2021. Sci. Total Environ.; 2023; 854, 158599. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2022.158599]

11. An, J.; Shi, Y.; Wang, J.; Zhu, B. Temporal variations of O₃ and NO_x in the urban background atmosphere of Nanjing, East China. Arch. Environ. Contam. Toxicol.; 2016; 71, pp. 224-234. [DOI: https://dx.doi.org/10.1007/s00244-016-0290-8]

12. Wang, W.-N.; Cheng, T.-H.; Gu, X.-F.; Chen, H.; Guo, H.; Wang, Y.; Bao, F.-W.; Shi, S.-Y.; Xu, B.-R.; Zuo, X. Assessing spatial and temporal patterns of observed ground-level ozone in China. Sci. Rep.; 2017; 7, 3651. [DOI: https://dx.doi.org/10.1038/s41598-017-03929-w]

13. Zhang, A.; Lin, J.; Chen, W.; Lin, M.; Lei, C. Spatial-Temporal Distribution Variation of Ground-Level Ozone in China’s Pearl River Delta Metropolitan Region. Int. J. Environ. Res. Public. Health; 2021; 18, 872. [DOI: https://dx.doi.org/10.3390/ijerph18030872] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33498400]

14. Tang, G.; Wang, Y.; Li, X.; Ji, D.; Gao, X. Spatial-temporal variations of surface ozone and ozone control strategy for Northern China. Atmos. Chem. Phys. Discuss.; 2011; 11, pp. 26057-26109.

15. Gao, C.; Xiu, A.; Zhang, X.; Chen, W.; Liu, Y.; Zhao, H.; Zhang, S. Spatiotemporal characteristics of ozone pollution and policy implications in Northeast China. Atmos. Pollut. Res.; 2020; 11, pp. 357-369. [DOI: https://dx.doi.org/10.1016/j.apr.2019.11.008]

16. Liu, H.; Liu, J.; Liu, Y.; Yi, K.; Yang, H.; Xiang, S.; Ma, J.; Tao, S. Spatiotemporal variability and driving factors of ground-level summertime ozone pollution over eastern China. Atmos. Environ.; 2021; 265, 118686. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2021.118686]

17. Lyu, Y.; Ju, Q.; Lv, F.; Feng, J.; Pang, X.; Li, X. Spatiotemporal variations of air pollutants and ozone prediction using machine learning algorithms in the Beijing-Tianjin-Hebei region from 2014 to 2021. Environ. Pollut.; 2022; 306, 119420. [DOI: https://dx.doi.org/10.1016/j.envpol.2022.119420] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35526642]

18. Yang, G.; Liu, Y.; Li, X. Spatiotemporal distribution of ground-level ozone in China at a city level. Sci. Rep.; 2020; 10, 7229. [DOI: https://dx.doi.org/10.1038/s41598-020-64111-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32350319]

19. Xue, T.; Zheng, Y.; Geng, G.; Xiao, Q.; Meng, X.; Wang, M.; Li, X.; Wu, N.; Zhang, Q.; Zhu, T. Estimating spatiotemporal variation in ambient ozone exposure during 2013–2017 using a data-fusion model. Environ. Sci. Technol.; 2020; 54, pp. 14877-14888. [DOI: https://dx.doi.org/10.1021/acs.est.0c03098]

20. Li, Y.; Shi, G.; Chen, Z. Spatial and temporal distribution characteristics of ground-level nitrogen dioxide and ozone across China during 2015–2020. Environ. Res. Lett.; 2021; 16, 124031. [DOI: https://dx.doi.org/10.1088/1748-9326/ac3794]

21. Liu, R.; Ma, Z.; Liu, Y.; Shao, Y.; Zhao, W.; Bi, J. Spatiotemporal distributions of surface ozone levels in China from 2005 to 2017: A machine learning approach. Environ. Int.; 2020; 142, 105823. [DOI: https://dx.doi.org/10.1016/j.envint.2020.105823]

22. Wei, J.; Li, Z.; Li, K.; Dickerson, R.R.; Pinker, R.T.; Wang, J.; Liu, X.; Sun, L.; Xue, W.; Cribb, M. Full-coverage mapping and spatiotemporal variations of ground-level ozone (O₃) pollution from 2013 to 2020 across China. Remote Sens. Environ.; 2021; 270, 112775.

23. Chang, C.-Y.; Faust, E.; Hou, X.; Lee, P.; Kim, H.C.; Hedquist, B.C.; Liao, K.-J. Investigating ambient ozone formation regimes in neighboring cities of shale plays in the Northeast United States using photochemical modeling and satellite retrievals. Atmos. Environ.; 2016; 142, pp. 152-170. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2016.06.058]

24. Choi, Y.; Kim, H.; Tong, D.; Lee, P. Summertime weekly cycles of observed and modeled NO_x and O₃ concentrations as a function of satellite-derived ozone production sensitivity and land use types over the Continental United States. Atmos. Chem. Phys.; 2012; 12, pp. 6291-6307. [DOI: https://dx.doi.org/10.5194/acp-12-6291-2012]

25. Duncan, B.N.; Yoshida, Y.; Olson, J.R.; Sillman, S.; Martin, R.V.; Lamsal, L.; Hu, Y.; Pickering, K.E.; Retscher, C.; Allen, D.J. Application of OMI observations to a space-based indicator of NO_x and VOC controls on surface ozone formation. Atmos. Environ.; 2010; 44, pp. 2213-2223. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2010.03.010]

26. Jin, X.; Fiore, A.; Boersma, K.F.; Smedt, I.D.; Valin, L. Inferring Changes in Summertime Surface Ozone–NO_x–VOC Chemistry over US Urban Areas from Two Decades of Satellite and Ground-Based Observations. Environ. Sci. Technol.; 2020; 54, pp. 6518-6529. [DOI: https://dx.doi.org/10.1021/acs.est.9b07785]

27. Pusede, S.; Cohen, R. On the observed response of ozone to NO_x and VOC reactivity reductions in San Joaquin Valley California 1995–present. Atmos. Chem. Phys.; 2012; 12, pp. 8323-8339.

28. Chi, X.; Liu, C.; Xie, Z.; Fan, G.; Wang, Y.; He, P.; Fan, S.; Hong, Q.; Wang, Z.; Yu, X. Observations of ozone vertical profiles and corresponding precursors in the low troposphere in Beijing, China. Atmos. Res.; 2018; 213, pp. 224-235. [DOI: https://dx.doi.org/10.1016/j.atmosres.2018.06.012]

29. Li, L.; Xie, S.; Zeng, L.; Wu, R.; Li, J. Characteristics of volatile organic compounds and their role in ground-level ozone formation in the Beijing-Tianjin-Hebei region, China. Atmos. Environ.; 2015; 113, pp. 247-254. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2015.05.021]

30. Wei, W.; Lv, Z.F.; Li, Y.; Wang, L.T.; Cheng, S.; Liu, H. A WRF-Chem model study of the impact of VOCs emission of a huge petro-chemical industrial zone on the summertime ozone in Beijing, China. Atmos. Environ.; 2018; 175, pp. 44-53. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2017.11.058]

31. An, J.; Zou, J.; Wang, J.; Lin, X.; Zhu, B. Differences in ozone photochemical characteristics between the megacity Nanjing and its suburban surroundings, Yangtze River Delta, China. Environ. Sci. Pollut. Res.; 2015; 22, pp. 19607-19617. [DOI: https://dx.doi.org/10.1007/s11356-015-5177-0]

32. Geng, F.; Tie, X.; Xu, J.; Zhou, G.; Peng, L.; Gao, W.; Tang, X.; Zhao, C. Characterizations of ozone, NO_x, and VOCs measured in Shanghai, China. Atmos. Environ.; 2008; 42, pp. 6873-6883. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2008.05.045]

33. Lu, H.; Lyu, X.; Cheng, H.; Ling, Z.; Guo, H. Overview on the spatial-temporal characteristics of the ozone formation regime in China. Environ. Sci. Process. Impacts; 2019; 21, pp. 916-929. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31089656]

34. Li, M.; Song, Y.; Mao, Z.; Liu, M.; Huang, X. Impacts of thermal circulations induced by urbanization on ozone formation in the Pearl River Delta region, China. Atmos. Environ.; 2016; 127, pp. 382-392. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2015.10.075]

35. Liu, P.; Song, H.; Wang, T.; Wang, F.; Li, X.; Miao, C.; Zhao, H. Effects of meteorological conditions and anthropogenic precursors on ground-level ozone concentrations in Chinese cities. Environ. Pollut.; 2020; 262, 114366. [DOI: https://dx.doi.org/10.1016/j.envpol.2020.114366] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32443214]

36. Lu, X.; Hong, J.; Zhang, L.; Cooper, O.R.; Schultz, M.G.; Xu, X.; Wang, T.; Gao, M.; Zhao, Y.; Zhang, Y. Severe surface ozone pollution in China: A global perspective. Environ. Sci. Technol. Lett.; 2018; 5, pp. 487-494. [DOI: https://dx.doi.org/10.1021/acs.estlett.8b00366]

37. Lu, X.; Zhang, L.; Chen, Y.; Zhou, M.; Zheng, B.; Li, K.; Liu, Y.; Lin, J.; Fu, T.-M.; Zhang, Q. Exploring 2016–2017 surface ozone pollution over China: Source contributions and meteorological influences. Atmos. Chem. Phys.; 2019; 19, pp. 8339-8361.

38. Wang, T.; Xue, L.; Brimblecombe, P.; Lam, Y.F.; Li, L.; Zhang, L. Ozone pollution in China: A review of concentrations, meteorological influences, chemical precursors, and effects. Sci. Total Environ.; 2017; 575, pp. 1582-1596.

39. Li, K.; Jacob, D.J.; Liao, H.; Shen, L.; Zhang, Q.; Bates, K.H. Anthropogenic drivers of 2013–2017 trends in summer surface ozone in China. Proc. Natl. Acad. Sci. USA; 2019; 116, pp. 422-427. [DOI: https://dx.doi.org/10.1073/pnas.1812168116]

40. Li, J.; Han, X.; Jin, M.; Zhang, X.; Wang, S. Globally analysing spatiotemporal trends of anthropogenic PM2.5 concentration and population’s PM2.5 exposure from 1998 to 2016. Environ. Int.; 2019; 128, pp. 46-62. [DOI: https://dx.doi.org/10.1016/j.envint.2019.04.026]

41. Li, G.; Haining, R.; Richardson, S.; Best, N. Space-time variability in burglary risk: A Bayesian spatio-temporal modelling approach. Spat. Stat.; 2014; 9, pp. 180-191. [DOI: https://dx.doi.org/10.1016/j.spasta.2014.03.006]

42. Malash, G.F.; El-Khaiary, M.I. Piecewise linear regression: A statistical method for the analysis of experimental adsorption data by the intraparticle-diffusion models. Chem. Eng. J.; 2010; 163, pp. 256-263. [DOI: https://dx.doi.org/10.1016/j.cej.2010.07.059]

43. Besag, J.; York, J.; Mollié, A.; Besag, J.; York, J.; Mollié, A. Bayesian image restoration, with two applications in spatial statistics. Ann. Inst. Stat. Math.; 1991; 43, pp. 1-20. [DOI: https://dx.doi.org/10.1007/BF00116466]

44. Lunn, D.J.; Thomas, A.; Best, N.; Spiegelhalter, D. WinBUGS—A Bayesian modelling framework: Concepts, structure, and extensibility. Stat. Comput.; 2000; 10, pp. 325-337. [DOI: https://dx.doi.org/10.1023/A:1008929526011]

45. Anselin, L. Spatial Econometrics: Methods and Models; Springer Science & Business Media: Berlin, Germany, 1988.

46. Brueckner, J.K. Strategic interaction among governments: An overview of empirical studies. Int. Reg. Sci. Rev.; 2003; 26, pp. 175-188. [DOI: https://dx.doi.org/10.1177/0160017602250974]

47. Elhorst, P.; Piras, G.; Arbia, G. Growth and Convergence in a Multiregional Model with Space-Time Dynamics. Geogr. Anal.; 2010; 42, pp. 338-355. [DOI: https://dx.doi.org/10.1111/j.1538-4632.2010.00796.x]

48. Hans, C. Bayesian lasso regression. Biometrika; 2009; 96, pp. 835-845. [DOI: https://dx.doi.org/10.1093/biomet/asp047]

49. Tibshirani, R. Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Ser. B Methodol.; 1996; 58, pp. 267-288. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1996.tb02080.x]

50. Rosenbaum, P.R.; Rubin, D.B. The central role of the propensity score in observational studies for causal effects. Biometrika; 1983; 70, pp. 41-55. [DOI: https://dx.doi.org/10.1093/biomet/70.1.41]

51. Ding, D.; Xing, J.; Wang, S.; Chang, X.; Hao, J. Impacts of emissions and meteorological changes on China’s ozone pollution in the warm seasons of 2013 and 2017. Front. Environ. Sci. Eng.; 2019; 13, pp. 1-9.

52. Liu, J.; Wang, L.; Li, M.; Liao, Z.; Sun, Y.; Song, T.; Gao, W.; Wang, Y.; Li, Y.; Ji, D. Quantifying the impact of synoptic circulation patterns on ozone variability in northern China from April to October 2013–2017. Atmos. Chem. Phys.; 2019; 19, pp. 14477-14492. [DOI: https://dx.doi.org/10.5194/acp-19-14477-2019]

53. Yu, Y.; Wang, Z.; He, T.; Meng, X.; Xie, S.; Yu, H. Driving factors of the significant increase in surface ozone in the Yangtze River Delta, China, during 2013–2017. Atmos. Pollut. Res.; 2019; 10, pp. 1357-1364. [DOI: https://dx.doi.org/10.1016/j.apr.2019.03.010]

54. Wei, J.; Li, Z.; Chen, X.; Li, C.; Sun, Y.; Wang, J.; Lyapustin, A.; Brasseur, G.; Jiang, M.; Sun, L. et al. Separating daily 1 km PM2.5 inorganic chemical composition in China since 2000 via deep learning integrating ground, satellite, and model data. Environ. Sci. Technol.; 2023; [DOI: https://dx.doi.org/10.1021/acs.est.3c00272] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37114869]

55. Churkina, G.; Kuik, F.; Bonn, B.; Lauer, A.; Grote, R.d.; Tomiak, K.; Butler, T.M. Effect of VOC emissions from vegetation on air quality in Berlin during a heatwave. Environ. Sci. Technol.; 2017; 51, pp. 6120-6130. [DOI: https://dx.doi.org/10.1021/acs.est.6b06514] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28513175]

56. Haagen-Smit, A.J. Chemistry and physiology of Los Angeles smog. Ind. Eng. Chem.; 1952; 44, pp. 1342-1346. [DOI: https://dx.doi.org/10.1021/ie50510a045]

57. Ma, M.; Gao, Y.; Ding, A.; Su, H.; Liao, H.; Wang, S.; Wang, X.; Zhao, B.; Zhang, S.; Fu, P. Development and assessment of a high-resolution biogenic emission inventory from urban green spaces in China. Environ. Sci. Technol.; 2021; 56, pp. 175-184. [DOI: https://dx.doi.org/10.1021/acs.est.1c06170]