Spatiotemporal Variation in Ground Level Ozone and Its Driving Factors: A Comparative Study of Coastal and Inland Cities in Eastern China

Variations in marine and terrestrial geographical environments can cause considerable differences in meteorological conditions, economic features, and population density (PD) levels between coastal and inland cities, which in turn can affect the urban air quality. In this study, a five-year (2016–2020) dataset encompassing air monitoring (from the China National Environmental Monitoring Centre), socioeconomic statistical (from the Shandong Province Bureau of Statistics) and meteorological data (from the U.S. National Centers for Environmental Information, National Oceanic and Atmospheric Administration) was employed to investigate the spatiotemporal distribution characteristics and underlying drivers of urban ozone (O3) in Shandong Province, a region with both land and sea environments in eastern China. The main research methods included the multiscale geographically weighted regression (MGWR) model and wavelet analysis. From 2016 to 2019, the O3 concentration increased year by year in most cities, but in 2020, the O3 concentration in all cities decreased. O3 concentration exhibited obvious regional differences, with higher levels in inland areas and lower levels in eastern coastal areas. The MGWR analysis results indicated the relationship between PD, urbanization rate (UR), and O3 was greater in coastal cities than that in the inland cities. Furthermore, the wavelet coherence (WTC) analysis results indicated that the daily maximum temperature was the most important factor influencing the O3 concentration. Compared with NO, NO2, and NOx (NOx ≡ NO + NO2), the ratio of NO2/NO was more coherent with O3. In addition, the temperature, the wind speed, nitrogen oxides, and fine particulate matter (PM2.5) exerted a greater impact on O3 in coastal cities than that in inland cities. In summary, the effects of the various abovementioned factors on O3 differed between coastal cities and inland cities. The present study could provide a scientific basis for targeted O3 pollution control in coastal and inland cities.

Volatile organic compounds (VOCs) and nitrogen oxides (NO x ≡ NO + NO 2 ) are precursors of O 3 [27]. Many studies have explored the relationship, photochemical reaction mechanism, and influence factors between O 3 and its precursors [28][29][30][31][32][33]. Meteorological conditions are important factors of O 3 generation and transmission. Radiation enhancement, temperature increase, and sunshine time extension facilitate O 3 generation, while a high relative humidity reduces O 3 [3,[34][35][36]. In addition to meteorological conditions, there may also be a relationship between socioeconomic factors and the spatial distribution of air pollution [37]. Many studies have analyzed the correlation between PM 2.5 and socioeconomic factors [38,39]. However, there is little research on the socioeconomic factors of O 3 . Even more notably, few studies have examined the combined impact of the O 3 precursors, meteorological conditions, and socioeconomic factors on O 3 pollution from a spatiotemporal perspective. These influencing factors may have different effects in different regions, which is worth exploring.
The difference in thermal conditions between land and oceans results in distinct climates between inland and coastal areas. Coupled with differences in the geographical environment, there could occur a comprehensive effect leading to different economic development and population distributions between inland and coastal areas. Within the context of global warming, the contrast between land and ocean warming levels is increasing (land warming enhanced), which intensifies aerosol pollution [40]. Considering that the influence of various factors on the O 3 concentration in coastal and inland cities may differ, it is therefore necessary to evaluate these differences.
The multiscale geographically weighted regression (MGWR) model is an extension of the geographically weighted regression (GWR) model, which is commonly used in the analysis of spatial influencing factors [41]. The MGWR model considers that different spatial scales exert different effects of the action mode and intensity on the spatial relationship between the considered factors and dependent variables, so it can explore the spatial heterogeneity of influencing factors [41]. MGWR has been widely used in the environmental field, such as PM 2.5 , PM 10 , NO 2 , and other air pollutants [42][43][44]. Zhan et al. (2022) has examined the spatial heterogeneity effects of both socioeconomic factors and natural factors on continuous air pollution with a MGWR [45]. Thus, we attempted to analyze the underlying drivers affecting the spatial distribution of O 3 and its spatial heterogeneity from socioeconomic aspects using MGWR.
Wavelet analysis was developed based on the Fourier transform [46]. The wavelet transform extends time series into the time-frequency space to overcome the limitations of the Fourier transform [47]. Wavelet analysis is a powerful analysis tool highly suitable to study nonstationary processes in the finite space-time domain [46,48]. This method has been widely applied in geophysics, economics, and public health research [49][50][51][52][53][54][55][56]. Therefore, this study will explore the potential drivers affecting the spatial distribution of O 3 from the perspective of meteorological conditions and other air pollutants through wavelet analysis (wavelet coherence (WTC)).
Shandong Province is located in eastern China, bordering the Bohai Sea and Yellow Sea. There occur both typical coastal cities and typical inland cities due to the large east-west span. In addition, Shandong Province ranks third in terms of its economy and second in terms of its population in China. The well-developed economy, large population, and high consumption of resources have caused serious air pollution, including O 3 pollution, in this area. Based on the above, we reasonably assume that the differences in the socioeconomic and natural environment of Shandong Province lead to significant differences in O 3 concentrations between coastal and inland areas, and the influence of various factors in different cities is different. This study will first attempt to combine MGWR and wavelet analysis to explore the influence of various factors on O 3 from multiple perspectives, so as to identify the key factors influencing the O 3 distribution. This study can provide the Chinese government with a more targeted reference for the treatment of O 3 pollution in coastal and inland cities.  (Figure 1). The region contains a complex, diverse terrain, and a long coastline, with four typical coastal cities: Qingdao, Yantai, Weihai and Rizhao. Shandong Province exhibits four distinct seasons, belonging to the warm temperate monsoon climate. The weather is changeable in spring, with less rain and windy sand-prone conditions. In summer, controlled by the southeast marine monsoon, southerly winds prevail. In autumn, the weather is sunny and moderate. Winter is controlled by the continental monsoon climate, mostly involving northerlies. The mountains in the central part of the territory are raised, which comprise mountains with an altitude higher than a kilometer, such as Mount Tai, Mount Lu, Mount Yi, and Mount Meng. The surrounding areas gradually transition from low mountains and hills into plains. Chinese government with a more targeted reference for the treatment of O3 pollution in coastal and inland cities.

Study Area
Shandong Province is located between 34°22.90′~38°24.01′ N and 114°47.50′~122° 42.30′ E (Figure 1). The region contains a complex, diverse terrain, and a long coastline, with four typical coastal cities: Qingdao, Yantai, Weihai and Rizhao. Shandong Province exhibits four distinct seasons, belonging to the warm temperate monsoon climate. The weather is changeable in spring, with less rain and windy sand-prone conditions. In summer, controlled by the southeast marine monsoon, southerly winds prevail. In autumn, the weather is sunny and moderate. Winter is controlled by the continental monsoon climate, mostly involving northerlies. The mountains in the central part of the territory are raised, which comprise mountains with an altitude higher than a kilometer, such as Mount Tai, Mount Lu, Mount Yi, and Mount Meng. The surrounding areas gradually transition from low mountains and hills into plains.

O3 and Other Air Pollutants Data
The DMA8 O3 concentration and the 24-h average concentrations of NO, NO2, NOx, PM2.5, and respirable particles (PM10) were considered in this study. Data pertaining to 16 cities in Shandong Province from 1 January 2016, to 31 December 2020, were obtained from the China National Environmental Monitoring Centre (CNEMC, http://www.cnemc.cn/, accessed on 13 September 2021). The annual O3 concentration is defined as the 90th percentile of DMA8 O3 concentration [57]. O3 monitoring and recording were carried out in strict accordance with Chinese national standards [58]. The data statistical validity of each evaluation project was implemented in accordance with the relevant provisions of the Chinese national standard [59].

O 3 and Other Air Pollutants Data
The DMA8 O 3 concentration and the 24-h average concentrations of NO, NO 2 , NO x , PM 2.5 , and respirable particles (PM 10 ) were considered in this study. Data pertaining to 16 cities in Shandong Province from 1 January 2016, to 31 December 2020, were obtained from the China National Environmental Monitoring Centre (CNEMC, http://www.cnemc.cn/, accessed on 13 September 2021). The annual O 3 concentration is defined as the 90th percentile of DMA8 O 3 concentration [57]. O 3 monitoring and recording were carried out in strict accordance with Chinese national standards [58]. The data statistical validity of each evaluation project was implemented in accordance with the relevant provisions of the Chinese national standard [59].

Socioeconomic and Meteorological Data
The socioeconomic data were derived from the statistics, including data on the population, gross domestic product (GDP), industrial power consumption, and number of civil vehicles (NCV) for 16 cities from 2016 to 2020 [60][61][62][63][64]. We also compiled the daily meteorological data from National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Information (https://www.ncei.noaa.gov/maps/global-summaries/, accessed on 5 September 2021). The meteorological monitoring sites are located in Jinan (116.9833 • E, 36.6833 • N) and Qingdao (120.3744 • E, 36.2661 • N). The specific monitoring items include average wind speed (WDSP, m/s), maximum continuous wind speed (MXSPD, m/s), maximum temperature (MAX, • C), and average temperature (TEMP, • C).

Methods and Mechanism
We used ArcGIS to visualize the spatiotemporal variation characteristics of O 3 concentration in Shandong Province from 2016 to 2020 and conduct hot spot analysis. OLS analysis was performed in R 4.0.5 invented by Rose Ihaka and Robert Gentleman of New Zealand. The GWR and MGWR models were implemented in MGWR 2.2 provided by the School of Geographical Sciences and Urban Planning at Arizona State University (https://sgsup.asu.edu accessed on 5 September 2021). Wavelet analysis was performed using publicly available MATLAB code.

GWR and MGWR
The GWR model is a new method to incorporate spatial correlation into a regression model, which allows regression coefficients to vary in space [41]. The expression of GWR model is as follows: where (u i , v i ) denotes the spatial coordinates of the ith observation point, β (u i , v i ) is the regression coefficient of the jth independent variable at the ith observation point, β 0 (u i , v i ) is the intercept of the model at the ith observation point, and e i is the error term. The MGWR model is an extension of the GWR model [41]. The largest difference between these models is the bandwidth. All variables in the GWR model exhibit the same bandwidth. However, the MGWR model specifies a dedicated bandwidth for each variable, thus reducing estimation errors and ensuring a more realistic and useful spatial process model. Different bandwidths can reveal the scale effect of different factors on O 3 concentration changes. Generally, the larger the bandwidth is, the lower the spatial heterogeneity [44]. The calculation equation is as follows: where (u i , v i ) denotes the coordinates of position i, bwj denotes the bandwidth considered by the regression coefficient of the jth variable, and β bw j (u i , v i ) is the regression coefficient of the jth variable at i.

Wavelet Analysis
Wavelet analysis is based on Fourier analysis. Elucidation of the localization characteristics of the analyzed object in the time and frequency domains constitutes the advantage of wavelet analysis [65]. The continuous wavelet transform (CWT) method entails wavelet superposition of different scales and displacement levels [66]. One frequently applied wavelet function is the Morlet wavelet, which is also adopted in this paper. Based on two CWTs, Grinsted et al. (2004) constructed the cross-wavelet transform (XWT), which exposes the high common power and relative phase in time-frequency space [67]. The XWT of two time series x n and y n is denoted as W XY , and the cross-wavelet power is denoted as |W XY |.
The complex argument arg(W XY ) can be interpreted as the local relative phase between x n and y n in time-frequency space and is given as follows: is the matrix of the smoothed cross-wavelet power spectra between x n and y n and Im and Re denote the imaginary and real parts, respectively, of . The phase angle can be determined to analyze the variation relationship (lag or consistent change) between two time series [49]. Consequently, left-and right-pointing phase angles indicate antiphase and in-phase relationships, respectively.
WTC analysis can be employed to study the correlation between two data series in time-frequency space. The WTC can significantly enhance the linear relationship and determine the covariance intensity between two time series [49,51]. Grinsted et al. (2004) defined the WTC of two time series as follows: where S is a smoothing operator and W X and W Y are the CWTs of x n and y n , respectively. R 2 ranges from 0 to 1.
where S scale denotes the smoothing level along the wavelet scale axis and S time denotes the smoothing level in time.

Choice of Two Typical Cities
Qingdao is the most economically well-developed city in Shandong Province, and Jinan is the capital city of Shandong Province. Furthermore, Qingdao is a coastal city, while Jinan is an inland city ( Figure 1). There are great geographical differences between these cities. Thus, this study chose Jinan and Qingdao as typical representative cities to explore the differences of the time-frequency relationship between O 3 and meteorological factors and other air pollutants between coastal and inland cities.

Spatiotemporal Distribution Characteristics of O 3
The annual O 3 concentration changes in 16 cities in Shandong Province are shown in year. An important reason for the observed increase could be enhancement in anthropogenic precursors. However, in 2020, the O 3 concentration was reduced due to strict control measures implemented by the Chinese government and the COVID-19 pandemic. The same trend is also shown in Figure 3, which indicates the annual changes in the O 3 concentration attainment rate in each city. The standard II concentration limit of the DMA8 O 3 concentration is 100 µg/m 3 , which is consistent with the air quality guideline (AQG) level of the World Health Organization Global Air Quality Guidelines (WHOGAQG) [68]. The standard II concentration limit of the DMA8 O 3 concentration is 160 µg/m 3 , which is consistent with the interim target 1 level of the WHOGAQG [68]. From 2016 to 2020, the proportion of the number of days when the O 3 concentration reached standard I in Shandong Province was 60.55%, 55   Moreover, there were obvious differences in the O3 distribution between inland and coastal areas. The O3 concentration in the inland areas (central and northwestern Shandong) was higher than that on the Jiaodong Peninsula (mainly coastal cities), which is consistent with Yao et al. (2019) [69]. This study further explored the spatial distribution characteristics of O3 through hot spot analysis (Figure 4). O3 cold spot areas were located in Qingdao, Yantai and Weihai. In addition, O3 high risk areas expanded from Jinan in 2016 to Jinan and surrounding areas in 2020, including Dezhou, Binzhou and Zibo. standard Ⅰ and standard Ⅱ, respectively. Only three coastal cities (Qingdao: 207/365; Rizhao: 188/365; Yantai: 187/365) exhibited more than half of the total number of days with the concentration reaching standard Ⅰ.  Moreover, there were obvious differences in the O3 distribution between inland and coastal areas. The O3 concentration in the inland areas (central and northwestern Shandong) was higher than that on the Jiaodong Peninsula (mainly coastal cities), which is consistent with Yao et al. (2019) [69]. This study further explored the spatial distribution characteristics of O3 through hot spot analysis (Figure 4). O3 cold spot areas were located in Qingdao, Yantai and Weihai. In addition, O3 high risk areas expanded from Jinan in 2016 to Jinan and surrounding areas in 2020, including Dezhou, Binzhou and Zibo.

Socioeconomic Impacts
In this study, redundant factors were eliminated according to the variance inflation factor (VIF) value (VIF < 7.5) of the ordinary least squares (OLS) model. VIF is a multicollinearity test method. The closer the VIF value is to 1, the lighter the multicollinearity is, and vice versa. Seven factors were selected: time, population density (PD, 10,000 peo-ple·hm −2 ), urbanization rate (UR, %), proportion of the secondary industry (PSI, %), output value of farming, forestry, animal husbandry and fishery (OFFAF, CNY), industrial power consumption (IPC, Billion kW·h), and NCV (unit). We tried to use PD to describe the population distribution, UR to represent the level of regional economic development, PSI to show the industrial structure, OFFAF to represent the development level of the primary industry, IPC to mean the industrial emissions of ozone, and NCV to represent the motor vehicle emissions of ozone. Three models, namely, OLS, GWR, and MGWR models, were employed to explore the socioeconomic relationship between these factors and the O3 concentration.
According to the corrected AICc, R 2 and adj-R 2 values of these three models (Table 1), the MGWR model achieved the best fitting effect. Compared with the OLS model, the R 2 and adj-R 2 values of the GWR and MGWR were greatly improved. But MGWR can more truly reflect the relationship between socioeconomic factors and O3 than GWR because the bandwidth of the MGWR model is variable, while the bandwidth of the GWR model is fixed. The MGWR model can reflect the spatial difference of the relationship between the independent variable and the dependent variable.

Socioeconomic Impacts
In this study, redundant factors were eliminated according to the variance inflation factor (VIF) value (VIF < 7.5) of the ordinary least squares (OLS) model. VIF is a multicollinearity test method. The closer the VIF value is to 1, the lighter the multicollinearity is, and vice versa. Seven factors were selected: time, population density (PD, 10,000 people·hm −2 ), urbanization rate (UR, %), proportion of the secondary industry (PSI, %), output value of farming, forestry, animal husbandry and fishery (OFFAF, CNY), industrial power consumption (IPC, Billion kW·h), and NCV (unit). We tried to use PD to describe the population distribution, UR to represent the level of regional economic development, PSI to show the industrial structure, OFFAF to represent the development level of the primary industry, IPC to mean the industrial emissions of ozone, and NCV to represent the motor vehicle emissions of ozone. Three models, namely, OLS, GWR, and MGWR models, were employed to explore the socioeconomic relationship between these factors and the O 3 concentration.
According to the corrected AICc, R 2 and adj-R 2 values of these three models (Table 1), the MGWR model achieved the best fitting effect. Compared with the OLS model, the R 2 and adj-R 2 values of the GWR and MGWR were greatly improved. But MGWR can more truly reflect the relationship between socioeconomic factors and O 3 than GWR because the bandwidth of the MGWR model is variable, while the bandwidth of the GWR model is fixed. The MGWR model can reflect the spatial difference of the relationship between the independent variable and the dependent variable. According to the OLS results (Table 1), there was a significant correlation between Time, PD, UR, PSI, and O 3 . Combined with Figures 2 and 3, the O 3 concentration had a distinct increasing trend from 2016 to 2020. The effects of PD and PSI on the O 3 concentration were positive across the whole province. An increase in PD and PSI could lead to an increase in the O 3 concentration to varying degrees. Densely populated areas contributed more to O 3 due to the consumption of many resources as fuel [70]. Cities with high industrial proportions produced more O 3 precursors, resulting in higher O 3 concentration. In most areas of Shandong Province, especially in coastal areas, there existed a stronger negative correlation between UR and the O 3 concentration. This is because coastal cities had a high degree of urbanization, but their O 3 concentration was low because of low emissions and favorable diffusion conditions. OFFAF, IPC, and NCV slightly affected O 3 and the relationship between them was not significant. Figure 5 indicates that compared with vehicle exhaust emissions, industrial emissions contributed more in the south of Shandong Province. The socioeconomic environment had a certain and extremely complex impact on the formation of O 3 . Greater population density and higher industrialization levels may lead to more O 3 precursors. In addition, due to the implementation of air pollution prevention and control policies, the areas with higher economic development levels are willing to invest higher more in air pollution control and better treatment effects. Notably, Lee et al. (2021) found that the number of houses, number of parking lots, and area of reconstruction projects were positively correlated with O 3 rates [71].
In this study, the effects of the mode and intensity of Time, PD, PSI, and NCV on the change in O 3 concentration were roughly similar on a large spatial scale, and the spatial relationship tended to remain stable (bandwidth is 71). The spatial heterogeneity in the impact of UR, OFFAF, and IPC on the O 3 concentration was high (bandwidths are 62, 57, and 57 respectively). What is more, Figure 5 shows the relationship between PD, UR, and O 3 is greater in coastal cities than that in the inland cities.
As a secondary air pollutant, processes leading to variations in surface O 3 are complex because they are impacted by natural and human factors. The impact of socioeconomic conditions on O 3 is only a small part, while the meteorological environment has a greater impact on O 3 . Therefore, it is difficult to explain the O 3 variations simply by socioeconomic statistics. To this end, we used meteorological and other air pollution data in the next section to further analyze its impact on O 3 .

Wavelet Power Spectrum of O3 in Jinan and Qingdao
The area enclosed by the black line in the Figure 6 indicates that the periodicity within the time series is significant. There occurred an obvious annual cycle in Qingdao and Jinan. The O3 concentrations in those two cities frequently fluctuated in summer (June to August) and autumn (September to November). In addition, the Qingdao O3 concentration exhibited a moderate cycle of approximately 150 days from 2016 to 2019. This indicates that double peaks occurred in the O3 concentration change in Qingdao, one peak in May or June and the other peak in September (Supplementary Figure S1). Affected by the summer monsoon climate, Qingdao is rainy in July and August. This led to a significant decline in the O3 concentration in summer due to lack of light [34,72].

Wavelet Power Spectrum of O 3 in Jinan and Qingdao
The area enclosed by the black line in the Figure 6 indicates that the periodicity within the time series is significant. There occurred an obvious annual cycle in Qingdao and Jinan. The O 3 concentrations in those two cities frequently fluctuated in summer (June to August) and autumn (September to November). In addition, the Qingdao O 3 concentration exhibited a moderate cycle of approximately 150 days from 2016 to 2019. This indicates that double peaks occurred in the O 3 concentration change in Qingdao, one peak in May or June and the other peak in September (Supplementary Figure S1). Affected by the summer monsoon climate, Qingdao is rainy in July and August. This led to a significant decline in the O 3 concentration in summer due to lack of light [34,72].

Wavelet Coherence Coefficient
In this study, NO, NO 2 , NO x , NO 2 /NO, PM 10 , PM 2.5 , WDSP, MXSPD, TEMP, and MAX were selected to explore the corresponding time-frequency relationship with O 3 . The arithmetic mean value of the WTC coefficient (ARsq) is provided in Table 2. In summary, the variation in O 3 in Qingdao was more closely related to other air pollutants than that in Jinan. Perhaps it is because Qingdao is greatly affected by the Yellow Sea, unlike Jinan, which is affected by many surrounding cities (Supplementary Figure S2). The time-frequency relationship between O 3 and other gas pollutants in Jinan was more complex. Among the meteorological elements, the coherence between the temperature and O 3 in Jinan was higher than that in Qingdao, while the coherence between the wind speed and O 3 in Qingdao was relatively higher. Previous studies have reported that O 3 generation largely depends on high temperatures and strong solar radiation [3,73]. Jinan's temperature fluctuated more greatly than that in Qingdao because Jinan is located far from the ocean. As a result, O 3 in Jinan was more vulnerable to the temperature. The wind speed notably impacted O 3 in Qingdao due to the prevalence of sea and land breezes. The daily maximum temperature exerted the strongest impact on O 3 in Jinan and Qingdao, followed by NO 2 /NO.  also believes that high temperature is the main meteorological driver in summer. In contrast to Jinan, PM 2.5 had a greater impact on O 3 in Qingdao [73].

Wavelet Coherence Coefficient
In this study, NO, NO2, NOx, NO2/NO, PM10, PM2.5, WDSP, MXSPD, TEMP, and MAX were selected to explore the corresponding time-frequency relationship with O3. The arithmetic mean value of the WTC coefficient (ARsq) is provided in Table 2. In summary, the variation in O3 in Qingdao was more closely related to other air pollutants than that in Jinan. Perhaps it is because Qingdao is greatly affected by the Yellow Sea, unlike Jinan, which is affected by many surrounding cities (Supplementary Figure S2). The time-frequency relationship between O3 and other gas pollutants in Jinan was more complex. Among the meteorological elements, the coherence between the temperature and O3 in Jinan was higher than that in Qingdao, while the coherence between the wind speed and O3 in Qingdao was relatively higher. Previous studies have reported that O3 generation largely depends on high temperatures and strong solar radiation [3,73]. Jinan's temperature fluctuated more greatly than that in Qingdao because Jinan is located far from the ocean. As a result, O3 in Jinan was more vulnerable to the temperature. The wind speed notably impacted O3 in Qingdao due to the prevalence of sea and land breezes. The daily maximum temperature exerted the strongest impact on O3 in Jinan and Qingdao, followed by NO2/NO.  also believes that high temperature is the main meteorological driver in summer. In contrast to Jinan, PM2.5 had a greater impact on O3 in Qingdao [73].

Time Frequency and Phase Relationship between O 3 and Other Air Pollutants
(1) NO, NO 2 , NO x , and NO 2 /NO. In the annual cycle, nitrogen oxides (including NO, NO 2 , and NO x ) and O 3 attained an inverse phase relationship, in which the value/phase of one parameter increased and that of the other parameter decreased (Figure 7). NO and O 3 exhibited an inverse phase relationship with a short period (less than 14 days). This may occur as a result of NO titration reaction (R1). NO 2 and O 3 , NO x and O 3 mostly attained a positive phase relationship in summer and an inverse phase relationship in winter in the short period. O 3 production in winter was generally in the NO x -saturated regime that high NO x concentration restrains O 3 formation. Simultaneously, high NO emissions titration contributed to a reduction in O 3 [11]. In summer, O 3 was in the NO x -limited regime that the increase of NO 2 was conducive to O 3 formation (R2). Xia et al. (2021) similarly showed that the daily O 3 concentration is higher when the daily NO 2 concentration is higher in summer [36]. In addition to NO x , VOCs are also important precursors for O 3 formation, which readily affect the variation of O 3 concentration in winter [74][75][76]. However, VOCs were not regarded as an influencing factor in this paper. The unexplained part of the figure may be mainly attributed to VOCs. In order to further clarify the relationship between NO x and O 3 , wavelet coherence analysis of NO 2 /NO and O 3 was performed ( Figure 7). As can be seen from the Table 2 and Figure 7, the coherence between NO 2 /NO and O 3 was stronger than that of NO, NO 2 , and NO x . There was an obvious annual cycle between the NO 2 /NO ratio and O 3 . In the annual cycle, the NO 2 /NO ratio and O 3 concentration in Jinan changed at the same time, while the arrow in Qingdao gradually tilted downward, indicating that the NO 2 /NO ratio had fluctuated before O 3 in Qingdao since 2018. In the short cycle, NO 2 /NO ratio had strong coherence with O 3 from May to July, and both increased simultaneously. This is consistent with the findings of Shao et al. (2009) and Wang et al. (2019) that the O 3 concentration increases with the NO 2 /NO ratio [28,74]. In spring, O 3 shifted from the NO x -saturated regime to the NO x -limited regime due to NO x concentration decrease (NO 2 gradually decreased, while NO plummeted). The magnitude of each change of NO 2 concentration in spring and summer was greater than that of NO (Supplementary Figure S3). Therefore, the increase of NO 2 concentration not only promoted the formation of O 3 , but also led to the increase of NO 2 /NO ratio. Wang et al. (2018) believe that NO titration is generally more significant on high-O 3 days, which further leads to an increase in the NO 2 /NO ratio [77]. Similarly, the reduction of NO 2 concentration inhibits the formation of O 3 . During the transitional stage between the NO x -saturated regime and the NO x -limited regime (spring and autumn), the changes in NO 2 /NO ratio and O 3 were reversed. This phenomenon was prominent in Jinan during the transition from NO x -saturated to NO x -limited (April and May), while Qingdao was more significant in the transition from NO x -limited to NO xsaturated (October). The concentration of NO 2 decreased gradually in spring, resulting in the decrease of the NO 2 /NO ratio. O 3 was weakened by NO x limitation, and its concentration gradually increased under favorable meteorological conditions (temperature rise and light enhancement). Li et al. (2021) also reached the same conclusion that the increase in O 3 observed in spring is driven by the reduction of NO x emissions [11]. In October, the ozone concentration in Qingdao was still at a high level. Under the condition of high concentration of NO 2 , high temperature and strong radiation, the concentration of O 3 increases. The consumption of NO 2 reduces the NO 2 /NO ratio, so that NO 2 /NO showed the opposite relationship with O 3 .
(2) PM 10 and PM 2.5 . Particulate matter (PM) (including PM 10 and PM 2.5 ) fluctuate frequently in winter and spring in Jinan and Qingdao. The phase trend of PM with O 3 was inclined to the left during the annual cycle, indicating that PM and O 3 exhibited an inverse phase relationship, and the change in O 3 preceded that in PM (Figure 7). In the short term (1~30 days), PM 2.5 and O 3 revealed opposite phases in winter and spring but the same phase in summer. In winter and spring, the PM 2.5 concentration was high due to the increase in energy consumption (heating). An ultrahigh PM 2.5 concentration could lead to light radiation weakening, which could reduce O 3 production. In addition, PM 2.5 scavenges HO 2 and NO x radicals, resulting in O 3 reduction [78]. O 3 production is high in summer due to the high temperatures and strong solar irradiation. The PM concentration was low in summer. When the PM 2.5 concentration increases to a certain extent, this could enhance the photochemical reactions of O 3 [36]. PM 10 [3]. This may occur because the O 3 concentration was the highest in the afternoon, and the highest temperature throughout the day also occurred in the afternoon [3.36]. Compared with Jinan, Qingdao attained a lower coherence between O 3 and temperature.

Limitations
However, the present study also suffers certain uncertainties and limitations. It is widely acknowledged that VOCs are some of the most important precursors for O 3 generation. Unfortunately, due to the lack of available high-resolution VOC data, this study did not comprehensively evaluate the possible impact of VOCs on the spatiotemporal distribution of O 3 . The analyzed O 3 data encompassed the average values of urban internal monitoring data according to Chinese national standard HJ 663-2013, while the obtained meteorological data included single-point data. Therefore, the spatial matching results of the meteorological and O 3 data were not completely consistent, which could lead to deviations in the analysis results in regard to the time-frequency relationship between O 3 and meteorological parameters.  The thick contour enclosed regions of greater than 95% confidence for a red noise process, the cone of influence where edge effects might distort the picture is shown as a lighter shade. The legend on the right shows the intensity of the wavelet power. The relative phase relationship is shown as arrows (with in-phase pointing right, anti-phase pointing left, and O 3 leading factors by 90 • pointing straight down) [67].

Conclusions
The combination of MGWR and wavelet analysis is excellent in exploring the effects of O 3 distribution from multiple perspectives. This is an innovative aspect of this study, which has not been seen in previous studies. There were significant spatiotemporal differences in the O 3 distribution in Shandong Province. On the one hand, the O 3 concentration increased annually from 2016 to 2019 but declined in 2020 under the influence of China's strict O 3 pollution cleanup policies and the global COVID-19 pandemic. On the other hand, the O 3 concentration in inland areas (namely, central, and northwestern Shandong) was higher, whereas that in eastern coastal areas was lower. In general, PD, UR, the wind speed (WDSP and MXSDP), nitrogen oxides, and PM 2.5 generated a greater impact on O 3 in a coastal city (Qingdao, China) than that in an inland city (Jinan, China). Whether coastal or inland cities, the ratio of NO 2 /NO was more coherent with O 3 than that NO, NO 2 , and NO x .
In summary, the reasons for the observed increase in the O 3 concentration in winter include the following: 1. reduction in nitrogen oxides; 2. reduction in the particle concentration; 3. transport of exogenous pollutants (VOCs and O 3 ) under the action of wind; 4. temperature rise. The observed O 3 concentration increase in summer is attributed to the high temperature, PM refraction and increase in nitrogen oxides. And the reasons for the decrease in O 3 concentration in summer include rainfall, the NO titration reaction, and reduction of NO 2 .
Nevertheless, China is located in eastern Asia and the west coast of the Pacific Ocean, with many important coastal and inland cities. Therefore, clarification of the differences in the main driving factors of O 3 between coastal and inland cities could enhance our understanding of O 3 pollution and improve the performance of O 3 prediction models. In addition, it would enable the Chinese government better control O 3 pollution. Regional O 3 pollution is a comprehensive reflection of the environment, society, economy, and development, and strengthening the leadership role of the government is therefore of crucial importance.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijerph19159687/s1, Figure S1: Histogram of monthly average O 3 concentration changes from 2016 to 2020; Figure S2: ARsq among cities in Shandong province from 2016 to 2020. * is indicates passing the 95% significance level test; Figure

Conflicts of Interest:
The authors declare that they have no known competing financial interest or personal relationships that could have appeared to influence the work reported in this paper.