Global Moran’s I is the best known and used method to reflect the similarity of spatial adjacent or adjacent regional cell property values (Chen, 2021). Moran’s I is calculated as the following Eq. (1):

where n is the number of monitoring cities; x_i and x_j refer to the attribute values of city i and j, respectively; w_ij is the spatial weight matrix between the regional units i and j.

Z_I is used to test the significance of global Moran’s I, and the calculation formula is as follows:

where Z_I is the Z test value for the global Moran’s I; E_I and Var_I are the mathematical expectation and covariance of the global Moran’s I, respectively.

The global spatial autocorrelation index can only reflect the overall process or trend, while it cannot reveal local differences. Therefore, it cannot specifically reflect the correlation and correlation between a city and its neighboring city (Zhou et al., 2020). To more accurately grasp the aggregation and differentiation characteristics of O₃ spatial agglomeration, the local Moran’s I (Anselin, 2010) is presented by Eq. (3):

where S is the standard deviation; If I_i is positive and significant, the position i is hot spot (high value agglomeration), if I_i is negative and significant, the position i is cold point (low value agglomeration).

To more scientifically and reasonably reflect the population exposure risk of O₃ in the study area; this paper calculates the population-weighted O₃ concentration value of a single grid by using a grid calculator (Fu and Kan, 2004). The calculation formula is as follows:

where PWEL_i is the population-weighted O₃ concentration average value, μg·m⁻³; i is the number of grid cells; P_i is the number of population in the grid; C_i is the O₃ concentration in this grid/(μg·m⁻³).

The Pearson correlation coefficient test was used to determine the correlation between the O₃ concentration and the meteorological factors (Dong et al., 2021). Correlation coefficient r can be calculated using Eq. (5):

where X_i, Y_i are the values of two random variables X and Y with linear relations, respectively; i is the number of samples, i = 1, 2, …, n; S_X and S_Y are the standard differences of the variable X, Y, respectively.

To better quantify the contribution of each socioeconomic factor to the O₃ concentration variety, the collinearity between explanatory variables should be first eliminated (Guo et al., 2016). The variance inflation factor (VIF) test is a classical method used to test the probable multicollinearity (Zhao et al., 2016; Che et al., 2019). In this study, multicollinearity was judged by statistics (T), robust probability (P), and variance inflation factor (VIF). Judgment methods were as follows: the larger the T is, the more significant the representation is; the smaller the p value is, the more useful variable is; if 0 < VIF < 10, there is no multicollinearity, and vice versa.

When establishing econometric models with cross-section data, the impact of explanatory variables on the interpreted variables may be different between regions due to the complexity, autocorrelation, and variability that this data exhibits. The GWR model can address the problem, which assumes that the economic behavior between regions is spatially heterogeneous and more realistic (Wang S. J. et al., 2020). The GWR model is as follows:

where x, y are independent and dependent variables, respectively; k is the number of independent variables; j is sample point; ε is regression residue; β₀ (u_j, v_j) is the intercept; and β_i (u_j, v_j) is the regression coefficient, changing with the sample point location. Each local β_i (u_j, v_j) is used to estimate its adjacent spatial observations.

Figure 1 shows the national distribution and the daily average exceedance rate of urban MDA8-90% in 2015 and 2018. In 2015 and 2018, the average annual concentration range were 62–202 and 74–215 μg·m⁻³, with a total of 60 and 111 cities exceeding the standard value (160 μg·m⁻³), accounting for 17.8% (60/337) and 32.9% (111/337), respectively. O₃ pollution distribution showed similarity to a certain extent in 2015 and 2018 and heavily polluted area was mainly concentrated in eastern China, such as Liaoning, Shandong, Hebei, Henan, and Jiangsu. Specifically, compared with 2015, O₃ distribution was more concentrated and serious in 2018. The polluted cities in severely polluted areas increased and the polluted areas gradually spread, for example, Shanxi, Shaanxi, and Anhui also began to suffer O₃ pollution.

As shown in the Figure 1, the areas with high O₃ exceedance rate were mainly concentrated in the eastern coastal areas of China, and it gradually spread to the central region. In 2015, the exceedance rate range of O₃ concentration was 0.00–24.93%. Among the studied 337 cities, the exceedance rates of 28 cities were above 15%, and those of four cities were above 20%, and these cities mainly located in Jiangsu, Shandong, Hebei, and Beijing. In 2018, the exceedance rate range of O₃ concentration was 0.00–30.14%, and the exceedance rates of 67 cities were above 15%, and those of 32 cities were above 20%. Generally, the Chinese O₃ concentrations and pollution area in 2018 have both increased compared with that in 2015; hence, it is of great significance to explore the spatial agglomeration and seasonal distribution characteristics.

Figure 2 shows the seasonal distribution of O₃ concentration in Chinese cities. Generally, the O₃ pollution in spring and summer was more serious, and the pollution scope were wider than that in fall and winter. Specifically, the high incidence areas of O₃ pollution in spring were mainly concentrated in central and eastern provinces, including Jiangsu, Shandong, Hebei, and Henan. In summer, pollution had spread to the western provinces, including Shaanxi, Gansu, Qinghai, Sichuan, and Inner Mongolia. The O₃ pollution areas decreased in fall compared with spring and summer, and these areas were mainly concentrated in Shandong, Jiangsu, and Anhui. In winter, O₃ pollution areas was mainly concentrated in the central and western provinces, including Sichuan, Qinghai, and Gansu.

The annual spatial agglomeration feature of national O₃ concentration is shown in Figure 3. The distribution areas of cold and hot spots were similar in 2015 and 2018 to a certain extent, which indicated that Chinese cities have formed a relatively stable and continuous pollution area. Hot spots were mainly distributed in eastern and central provinces, including Shanghai, Jiangsu, Anhui, Shandong, Beijing, Hebei, Shanxi, and Henan. Moreover, compared with 2015, Shaanxi, Inner Mongolia, Hubei, and Zhejiang have gradually become hot agglomeration areas in 2018. Cold points were mainly distributed in Guangxi, Guangdong, Hainan, Xinjiang, and Heilongjiang. The causes of the distributions are regional transport of O₃ and similar large-scale meteorological conditions. This phenomenon suggests that joint efforts among urban agglomerations are crucial to control O₃ pollution in the region, rather than just to control O₃ emissions in individual cities.

Spatial agglomeration characteristics of the O₃ concentration in spring, summer, fall, and winter in 2015 and 2018 are shown in Figure 4. In terms of season, in spring and summer, hot spots were mainly distributed in eastern, northern, and central provinces, including Inner Mongolia, Liaoning, Beijing, Hebei, Shanxi, Shandong, Henan, Anhui, Jiangsu, Shanghai, and Shaanxi. Therefore, the spring and summer were a critical period for O₃ pollution control in these provinces. Cold spots were mainly distributed in southwestern provinces, including Sichuan, Chongqing, Hunan, Jiangxi, Fujian, Yunnan, Guizhou, Guangxi, Guangdong, and Hainan. In fall, hot spots expanded in the southeast and were mainly concentrated in Hebei, Beijing, Shanxi, Shandong, Henan, Anhui, Jiangsu, Zhejiang, Jiangxi, Fujian, and Guangdong. Cold spots areas were gradually concentrated in the central provinces, including Sichuan, Chongqing, Guizhou, and Yunnan. Overall, O₃ pollution scale was large and these areas mainly concentrated in the central and eastern urban agglomerations, and showed distinct seasonal characteristics. Therefore, it is necessary to deeply strengthen the joint control measures between urban agglomerations in the central and eastern regions, and strengthen seasonal regulation, especially in spring and summer.

Table 2 shows that the correlation between O₃ concentration and temperature in different seasons. The annual O₃ concentration of the three urban agglomerations was significantly positively correlated with the temperature and the correlation is significant in the different seasons. In addition, seasonal differences were obvious. The correlation coefficient of the three urban agglomerations in spring and autumn was significantly higher than that in summer and winter. Among them, the CP-UA failed the confidence (bilateral) significance test of 0.01 in the winter in 2018 and the summer in 2015. One reason is possible that O₃ photochemical reaction is inefficient due to the influence of other factors, such as weak solar radiation in winter, rainfall and wind speed in summer (Chen X. P. et al., 2019).

Table 3 shows the correlation between O₃ concentration and relative humidity in different seasons. The annual relative humidity of the three urban agglomerations were significantly negatively correlated with O₃ concentration. In the three urban agglomerations, the absolute value of the correlation coefficient increased from north to south. The correlation was the most significant in summer, and the ranking of absolute value of correlation coefficient increased from north to south, which was the same as the annual ranking. It was due to that when the atmospheric relative humidity increases, it is accompanied by an increase in cloud cover, which leads to an increase in precipitation. These meteorological conditions are not conducive to the formation and accumulation of O₃, which leads to the decrease of O₃ concentration (Liang et al., 2019; Bai et al., 2022). In addition, it may be that when the RH is high, the photochemical decomposition of water vapor will produce more reactive groups and react with O₃, which reduces the concentration of O₃ (Tan et al., 2007). And this inhibition effect was more significant in the summer of the central and northern regions. Among them, the seasonal correlation coefficient of BTH-UA fluctuated significantly. The CP-UA was summer > spring and fall > winter, and the correlation coefficient was quite different, especially in summer and winter. The YRD-UA was summer > fall > spring > winter in 2015, and fall > spring > summer > winter in 2018, and the correlation coefficient fluctuation was small.

Wind speed, especially near-ground wind speed, determines the speed of pollutant handling and dilution (Yang et al., 2021). The correlation between O₃ concentration and wind speed in different seasons was shown in Table 4. The annual O₃ concentration of the CP-UA and the YRD-UA showed significantly negative correlation with the wind speed. While annual O₃ concentration of the BTH-UA was significantly positively correlated with the wind speed. From the seasonal point of view, the O₃ concentration in spring, summer, and autumn of the three urban agglomerations all showed significantly negative correlation with wind speed, while the O₃ concentration in winter showed significantly positive correlation with wind speed. It might have the reasons of O₃ increasing were the elevation of atmospheric boundary height and the increase of vertical momentum transport due to the increase of wind speed, which then promotes the transfer of O₃ to the ground (Chen et al., 2022; Liu et al., 2022). The reason for the O₃ decrease may be that the wind speed increases the horizontal diffusion movement of O₃.

The meteorological factors alter the O₃ concentration through physical and chemical processes. Among the meteorological factors, temperature was most associated with O₃ concentration in 2015 and 2018 (Tables 2–4). Temperature was directly affecting O₃ concentration by affecting the photochemical reaction generation efficiency of the O₃ (Chen et al., 2017), higher temperature, more frequent molecular collisions, and accelerated photochemical reaction rates. The wind speed and relative humidity had a negative effect on the O₃ concentration in BTH-UA, YRD-UA, and CP-UA. It was due to that high wind speed can promote horizontal diffusion and reduce the accumulation of air pollutants in headwind areas (Wang et al., 2019), indicating that the enhanced atmospheric diffusions efficiently reduced O₃ concentration levels. Moreover, O₃ pollution has changed significantly along with the season, so the seasonal O₃ pollution response measures should be strengthened. First, continue to strengthen the O₃ response in summer, and spring and fall should also be valued. Specifically, avoid or reduce VOCs process production in March–November, and stagger the peak emission in the O₃-prone period (12:00–17:00). Second, formulate a positive inventory of seasonal VOCs intensive emission reduction measures and implemented differentiated emission reduction measures. Enterprises that conform to the conditions of the positive inventory may not implement strengthened emission reduction. Third, in the critical areas such as BTH-UA, YRD-UA, and CP-UA, formulate pollution control plans for key industries, such as painting, printing, and textile, to help enterprises effectively carry out comprehensive control of VOCs.

This study conducted multicollinearity tests of 14 socioeconomic factors of agglomeration regions in 2015 and 2018, including total population at year-end (TP), population density (PD), the gross domestic product (GDP), the per capita GDP (Per_GDP), share of the primary industry in GDP (%; Primary), share of the secondary industry in GDP (%; Secondary), share of the tertiary industry in the GDP (%; Tertiary), industrial soot emissions (ISE), annual average population (AP), total land area of administrative region (Administration), area of green land (Green), industrial SO₂ emission (ISO₂), average annual concentration of inhaled fine particulate matter (PM₁₀), and annual electricity consumption (AEC). The ArcGIS 10.2 analysis results were shown in Supplementary Table S2. The VIF of TP, GDP, Primary, Secondary, Tertiary, AP, and AEC were more than 10, which indicated that these seven factors had strong multicollinearity by themselves or with other factors, so these seven factors should be discarded. While PD, Per_GDP, ISE, Administration, Green, ISO₂, and PM₁₀ do not exist seriously multicollinearity, which can be used for GWR.

The GWR model was fitted with seven socioeconomic factors (PD, Per_GDP, ISE, Administration, Green, ISO₂, and PM₁₀). The number of conditions on the sample points in 2015 and 2018 were all lower than 30, showing that there was no local collinearity. The coefficients of variables were all significant at the level of 1%, and the model estimation results were credible. The R² and adjusted R² were 0.46 and 0.31 in 2015, 0.66 and 0.56 in 2018 respectively, indicating the model had a better fit for both metrics in 2018. The standardized residual distribution was shown in Figure 7. The residue in the regression residue spatial distribution diagram was spatially random, demonstrating that the regression residue obeyed to the normal distribution. The regression results of normalized residual ranges were −2.51 to 2.07 in 2015, −3.11 to 1.63 in 2018, and values in the range of −2.00 to 2.00 accounted for 92.42 and 93.94% of the total results respectively, revealing the GWR model fit well. To further analyze the spatial changes of the various socioeconomic factors and the significance level of the O₃, Supplementary Table S3 reports the descriptive statistical results of the various socioeconomic factors of the GWR regression model.

The impact of PD, Per_GDP, ISE, Administration, Green, ISO₂, and PM₁₀ on O₃ pollution of urban agglomeration could be measured by the estimated coefficient of the influencing factors in various regions derived from the GWR model. The greater the coefficient of factor is, the greater the influence is, and the positive and negative of the coefficient represents the directionality of this factor. The parameter of the seven variables selected by the GWR differed in each region, indicating a spatial variation in the influence of each variable on the O₃ concentration.

The urbanization factors included PD, Administration, and Green, and the adjusted R² were 0.42, 0.39, and 0.34 in 2015 respectively, and those were 0.57, 0.54, and 0.54 in 2018. This statistically demonstrated that the O₃ concentration distribution was closely related to the degree of urbanization. Figure 7 shows the distribution of the urbanization factors (PD, Administration, and Green) regression coefficient in 2015 and 2018.

Areas with high population density tend to have more pollution sources and pollution activity. Standardized residues for GWR between O₃ and PD ranged from −2.88 to 1.96 in 2015, and − 3.23 to 1.70 in 2018, values in the range of −2.00 to 2.00 accounted for 93.94 and 95.45% of the total results, respectively. Figure 7 shows that population density was positively associated with O₃ concentration in the majority cities, namely O₃ concentrations increased with the density of the population. The top five cities with the greatest correlation influence in 2015 were Chuzhou, Bengbu, Suzhou, Hefei, and Ma’anshan, located in the east of the CP-UA and the west of the YRD-UA. The top five cities with the highest correlation influence in 2018 were Qinhuangdao, Tangshan, Chengde, Xinyang, and Wuhu, mainly located in the BTH-UA. It was due to densely populated areas where human activity was more intense, and O₃ pollution was closely related to precursor emissions, such as VOCs, CO, and NOx. Therefore, the future control policies should be deeply concentrated in BTH-UA areas with high population density.

As shown in Figure 8, O₃ pollution was negatively correlated with Administration and Green in the majority of cities. Standardized residues for GWR between O₃ and Administration ranged from −2.85 to 2.12 in 2015, −3.72 to 2.32 in 2018, values in the range of −2.00 to 2.00 accounted for 92.42 and 92.42% of the total results, respectively. Standardized residues for GWR between O₃ and Green ranged from −3.05 to 1.70 in 2015, −2.84 to 1.74 in 2018, values in the range of −2.00 to 2.00 accounted for 93.94 and 95.45% of the total results, respectively. The high-value area was located in the YRD-UA, and the low-value area was located in the BTH-UA and the area between the border of the BTH-UA and CP-UA, which shows that the larger the green space area in this area is, the more it was conducive to reduce O₃ pollution. While some cities in the BTH-UA such as Tangshan, Qinhuangdao, Hengshui, Xingtai, Handan, and CP-UA such as Liaocheng, Puyang, Xinxiang, Zhengzhou, Administration, and Green were positively related to O₃ pollution in the majority cities. It was due to that the high O₃ concentration caused a certain degree of damage to plants, and greening did not play a role in reducing O₃ pollution. Overall, green space could help to improve air pollution and it was necessary to further expand green space in various cities.

To alleviate O3 pollution caused by urbanization, the government needs to enhance the role of spatial allocation in urban planning, and create a spatial intensive urban pattern. Although central cities with high population density help to give full play to the advantages of agglomeration economy, high population will increase the ecological and environmental pressure and drive O₃ pollution. The reasonable layout and design of urban architectural could promote the dispersions of O₃ and improve the air quality. With the acceleration of urbanization in China, there were a large number of construction activities in developed and developing cities. Therefore, reasonable development plan of the administrative area and green space will be conducive to reduce O₃ pollution and improve the urban air quality.

According to the analysis of GWR model, the adjusted R² of Per_GDP in 2015 and 2018 was 0.38 and 0.59, respectively (Supplementary Table S3). And as shown in Figure 9, O₃ pollution was positively correlated with Per_GDP in the majority cities. Standardized residues for GWR between O₃ and Per_GDP ranged from −2.76 to 1.66 in 2015, and − 2.83 to 1.71 in 2018, values in the range of −2.00 to 2.00 accounted for 93.94 and 95.45% of the total results, respectively. The high-value area was located in the YRD-UA, and the low-value area was located in the BTH-UA.

Yang (2021) argued that China’s economic growth and environmental pollution showed a “U-shaped” relationship, that was, when the economic level was low, pollution improves with economic growth, and when it reaches an “inflection point,” pollution will deteriorate with economic growth. In this study, with the rapid development of economy from 2015 to 2018, O₃ pollution deteriorated. It could be seen that the BTH-UA, CP-UA, and YRD-UA were all in this “U-shaped” climbing period. With the further growth of the economy, regional O₃ pollution had a trend to deteriorate. Per_GDP was an important indicator reflecting the level of regional economic development. To realize the sustainable development of environmental protection and economic growth, the BTH-UA, CP-UA, and YRD-UA (especially YRD-UA) needed to build a joint pollution prevention and control model dominated by economic coordination and supplemented by policy and management coordination.

In this study, the anthropogenic factors included ISE, ISO₂, and PM₁₀. The adjusted R² of ISE, ISO₂, and PM₁₀ were 0.34, 0.45, and 0.31 in 2015, and those were 0.56, 0.59, and 0.6 in 2018, respectively (Supplementary Table S3). Standardized residues for GWR between O₃ and ISE ranged from −3.25 to 1.73 and −2.81 to 1.57 in 2015, 2018, values in the range of −2.00 to 2.00 accounted for 95.45 and 93.94% of the total results, respectively. Standardized residues for GWR between O₃ and ISO₂ ranged from −2.80 to 1.79 in 2015, and − 2.83 to 1.78 in 2018, values in the range of −2.00 to 2.00 accounted for 93.94 and 96.97% of the total results, respectively. Standardized residues for GWR between O₃ and PM₁₀ ranged from −3.19 to 1.63, values in the range of −2.00 to 2.00 accounted for 93.94% in 2015, and ranged from −2.96 to 2.00 and values in the range of −2.00 to 2.00 accounted for 95.45% of the total results in 2018. O₃ pollution was positively correlated with the anthropogenic factors (ISE, ISO₂, and PM₁₀) in the majority cities as shown in Figure 10. It indicated that industrial soot emissions, industrial SO₂ emissions, and inhaled fine particulate matter increased O₃ pollution. The high-value area was located in the CP-UA and YRD-UA, and the low-value area was located in the BTH-UA.

Compared with 2015, the impacts of ISE on O₃ pollution in 2018 were reduced, while the impacts of ISO₂ and PM₁₀ emission were more significant. These mean that the upcoming industrial emission control policy should place greater emphasis on limiting the SO₂ and PM₁₀ emissions in the CP-UA and YRD-UA. Specifically, SO₂ emissions should be strictly controlled in the YRD-UA. And to control PM₁₀ emissions, all kinds of open-air incineration should be strictly controlled and the government should strengthen the main responsibility of CP-UA governments for straw burning at all levels. In addition, the major sources of pollution in SO₂, soot, and PM₁₀ were industrial emissions. Therefore, the government needs to grasp the development direction of environmental protection technology of the “Made in China 2025” strategy, and encourage industrial enterprises to research, develop and introduce environmental protection technology, and accelerate the adjustment and optimization of the industrial structure in the BTH-UA, CP-UA, and YRD-UA. It will help to reduce the environmental burden. Specifically, control of exhaust pollution from mobile pollution sources should be strengthened and rigid industrial emission standards should be established, factories, which could not meet the standards, should be thus eliminated.