Xiating Chen
School of Mathematics and Statistics, Guilin University of Technology, Guilin, China.
DOI: 10.4236/jamp.2026.143056   PDF    HTML   XML   62 Downloads   312 Views   Citations

Abstract

Under the context of urban compound air pollution, the joint dependence structure and spatiotemporal characteristics between NO₂ and PM_2.5 have not yet been systematically quantified. To reveal their co-variation mechanisms, this study utilizes daily average data from 31 air quality monitoring stations covering urban and suburban areas of Beijing during 2023-2024. A staged analytical framework integrating overall modeling, spatial clustering, and seasonal stratification is constructed. First, a mixture copula model is employed to characterize the overall joint distribution of NO₂ and PM_2.5. Subsequently, statistical fingerprint clustering combined with seasonal analysis is applied to identify spatial and temporal heterogeneity in the dependence structure. The results indicate a moderately strong positive dependence between NO₂ and PM_2.5 across Beijing, with an asymmetric structure featuring upper-tail-dominant dependence structure. This suggests that the two pollutants are more likely to increase synchronously under high-pollution conditions, implying a joint amplification effect during extreme pollution episodes. Spatial clustering analysis reveals significant gradient differences in concentration levels among monitoring stations, particularly for NO₂, which exhibits significant spatial stratification. However, the form of the dependence structure remains relatively stable, indicating that spatial location primarily influences pollution intensity rather than coupling mechanisms. Seasonal variation substantially restructures the joint distribution. During the heating season, the two pollutants exhibit strong positive dependence and lower tail dependence, whereas in the non-heating season, the dependence is characterized by upper-tail-dominant feature. These findings suggest that seasonal factors exert a stronger regulatory effect on the pollutant dependence structure than spatial differences. This study reveals the spatiotemporal characteristics of the synergistic evolution of NO₂ and PM_2.5 in Beijing, providing support for joint risk assessment of extreme pollution events and seasonally differentiated coordinated control strategies.

Keywords

Dependence Modeling, Mixture Copula, Air Pollution, Clustering, Spatiotemporal Characteristics

Share and Cite:

1. Introduction

With the accelerating pace of global urbanization, urban air pollution has become one of the core environmental challenges constraining sustainable development and public health security [1] [2]. Among various atmospheric pollutants, fine particulate matter (PM_2.5) and nitrogen dioxide (NO₂) are considered representative urban pollutants due to their widespread emission sources, relatively high ambient concentrations, and substantial health risks. Many epidemiological studies have demonstrated that exposure to NO₂ and PM_2.5 is strongly associated with respiratory diseases, cardiovascular diseases, and premature mortality [3]-[5]. In severe pollution episodes, these two pollutants often exhibit synchronous increases, with joint exposure effects potentially exceeding those of individual pollutants. Therefore, revealing the dependence structure between NO₂ and PM_2.5 is of great significance for understanding the formation mechanisms of compound air pollution and conducting joint risk assessments.

In recent years, copula models have been widely applied in multivariate dependence modeling of air pollution variables. As a flexible multivariate dependence modeling tool, copula can precisely describe the joint structure between variables without assuming marginal distribution forms, making it particularly suitable for air pollutant data with significant asymmetry and tail dependence [6]-[9]. In practical applications, Archimedean copulas are frequently employed to describe different types of tail dependence structures [10]-[12]. Previous studies have applied copula models to joint analyses of different regions and pollution combinations. For example, Masseran et al. [13] combined copula models with other statistical models to analyze the dependence structure among multiple air pollutants in Malaysia. Zhang et al. [11] employed five copula models to study the multiscale correlations between PM_2.5 and O₃ in three major Chinese cities. Abbasi et al. [10] integrated time-series clustering with classical Archimedean copula models to explore the temporal variation and dependence between PM_2.5 and associated pollutants in Tehran. Sak et al. [14] assessed PM_2.5 pollution risk in multiple Chinese cities based on copula models. Chan et al. [15] constructed multivariate copula models to investigate dependence relationships of extreme air pollution events.

However, a single copula model is typically capable of capturing only one specific type of tail dependence. In real atmospheric environments, pollutant dependence structures are influenced jointly by emission sources, meteorological conditions, and regional transport processes, often resulting in composite structural characteristics. To overcome this limitation, recent studies have introduced linear combinations of copula models (mixture copula models), which more comprehensively fit complex dependence structures by weighting different copulas [16]-[19]. The joint distribution of air pollutants frequently exhibits non-normality, asymmetry, and tail dependence characteristics [4] [10] [20], which traditional linear correlation analysis cannot fully reveal. Therefore, this study adopts a mixture copula model composed of Gumbel, Clayton and Frank copulas to flexibly capture the potential multiple dependence patterns between NO₂ and PM_2.5.

Compared to existing studies, this paper extends the methodological framework hierarchically. An overall dependence baseline model is first established, followed by the integration of spatial clustering, seasonal stratification, and copula modeling to progressively reveal the spatiotemporal heterogeneity of the joint distribution between NO₂ and PM_2.5. Through this stepwise analytical framework, the paper aims to uncover the structural evolution patterns of their joint distribution in Beijing, thereby providing quantitative support for refined air quality management and differentiated prevention strategies. The proposed framework also offers methodological reference for multi-station, multi-pollutant joint risk assessment in other metropolitan regions.

The structure of this paper is as follows. Section 2 introduces the study area and data. Section 3 describes the research methods and model fundamentals. Section 4 presents the model fitting results and provides comprehensive analysis and discussion. Section 5 summarizes the main conclusions.

2. Study Area and Data

2.1. Study Area and Data Sources

Beijing is located in the northern part of China and is characterized by a temperate continental monsoon climate. Influenced by terrain and urban heat island effects, atmospheric diffusion and pollutant accumulation processes in this region exhibit substantial complexity, making it a typical case for investigating compound urban air pollution mechanisms. As the political, economic, and cultural center of China, Beijing has a dense population, frequent economic activities, and high energy consumption intensity. Consequently, modeling air pollutant behavior in Beijing is of great significance.

The data for this study are sourced from the Beijing Municipal Environmental Monitoring Center. Initially 35 air quality monitoring stations within Beijing were selected (station distribution shown in Figure 1). These stations cover urban and suburban functional zones, comprehensively reflecting the spatial distribution characteristics of air quality across the city. The research data consist of daily average concentrations of NO₂ and PM_2.5 at each station from January 1, 2023, to December 31, 2024 (units: μg/m³).

Figure 1. Distribution of air quality monitoring stations in Beijing.

2.2. Data Preprocessing

In practical air quality monitoring, data missing is inevitable due to equipment maintenance, power interruptions, or instrument failures. Existing studies have indicated that missing air quality data typically follows a random missing pattern [21] [22]. If not handled properly, missing data may cause bias in statistical modeling and parameter estimation. To ensure data integrity and reliability of analysis results, a hierarchical preprocessing strategy is adopted. If a monitoring station has consecutive missing data exceeding 10 days, it is considered to have insufficient data integrity and is entirely excluded. For a small number of discrete missing values, a sliding window interpolation method is used [23].

After data screening and preprocessing, four stations—Huairou New Town, Shunyi New Town, Mentougou Sanjiadian, and Changping Nanshao—are excluded due to continuous missing records exceeding 10 days, leaving 31 monitoring stations for subsequent analysis. Statistical results indicate that the overall missing rates for NO₂ and PM_2.5 data across valid stations are 0.26% and 0.29%, respectively, with missing patterns primarily consisting of short-term discrete gaps (1 - 4 consecutive days).

To address this pattern, a sliding window linear interpolation approach is employed to impute missing values. For periods of consecutive missing days $m$ ($1\le m\le 4$ ), a window centered on the missing interval is constructed by taking $\lceil m/2\rceil +1$ valid observations both before and after the missing segment. Linear interpolation is performed based on the valid data within the window. Prior to interpolation, extreme concentration values of NO₂ and PM_2.5 are identified using the Z-score method ($\left|Z\right|>3$ ) and preserve without participating in the interpolation calculations. In addition, interpolated values are constrained within the 5th–95th percentile range of the corresponding pollutant at each station to avoid distortion of extreme concentration distribution. This preprocessing procedure maintains sample size while improving data quality and the stability of model estimation [23].

To characterize the overall pollution level in Beijing, daily average concentrations from multiple stations are aggregated to construct a regional-scale dataset. Specifically, the arithmetic mean of NO₂ and PM_2.5 concentrations across all valid stations is calculated for each observation day, serving as a regional representative concentration.

3. Methodology

3.1. Copula Model

To characterize the nonlinear and asymmetric dependence structure between NO₂ and PM_2.5, a joint distribution modeling framework based on copula theory is adopted. Copula models serve as a powerful tool that links multivariate joint distributions with their marginal distributions, enabling the separation of dependence structure from marginal behavior and allowing independent modeling of each component. According to Sklar’s theorem [6] [24], for two random variables $X$ and $Y$ with joint distribution function $H\left(x,y\right)$ and marginal distribution functions $F_X\left(x\right)$ and $F_Y\left(y\right)$ , respectively, there exists a copula function $C:{\left[0,1\right]}^2\to \left[0,1\right]$ such that

$H\left(x,y\right)=C\left(F_X\left(x\right),F_Y\left(y\right)\right).$

Different copulas can capture distinct tail dependence structures [25]. This study selects three classic Archimedean copulas—Gumbel, Clayton and Frank—primarily based on their theoretical suitability and relevance to atmospheric physics. Gumbel copula captures upper tail dependence, reflecting the joint increase of pollutant concentrations during severe pollution episodes. Clayton copula characterizes lower tail dependence and is capable of modeling the synchronous decrease of pollutant concentrations under clean-air conditions. Frank copula describes symmetric dependence, reflecting moderate correlations under normal pollution levels. Together, these three copulas provide complementary representations of dependence across the full concentration range and have been widely used in air pollution dependence modeling. Although other copula families could also be considered, their tail dependence characteristics largely overlap with those of the above models and would substantially increase model complexity. Since the Gumbel-Clayton-Frank combination flexibly fits the dependence features of the pollutants in this study while maintaining model simplicity and interpretability, additional copulas are not included. Therefore, this study adopts the following mixture copula [6] [17]:

$C_{\text{mix}}\left(u,v\right)={\alpha }_1{C}_G\left(u,v;{\theta }_1\right)+{\alpha }_2{C}_C\left(u,v;{\theta }_2\right)+{\alpha }_3{C}_F\left(u,v;{\theta }_3\right),$ (1)

where $C_G,C_C,C_F$ denote Gumbel, Clayton and Frank copulas, respectively; ${\theta }_1,{\theta }_2,{\theta }_3$ are the intrinsic parameters of Gumbel, Clayton and Frank copulas, satisfying ${\theta }_1\in \left[1,\infty \right)$ , ${\theta }_2\in \left(0,\infty \right)$ , ${\theta }_3\ne 0$ ; ${\alpha }_1,{\alpha }_2,{\alpha }_3$ are the weight parameters corresponding to the copulas in the mixture model, satisfying ${\alpha }_i\in \left[0,1\right]$ for $i=1,2,3$ and ${\displaystyle {\sum }_{i=1}^3{\alpha }_i}=1$ . The three copulas are given by:

$C_G=\exp \left(-{\left({\left(-\log u\right)}^{{\theta }_1}+{\left(-\log v\right)}^{{\theta }_1}\right)}^{1/{\theta }_1}\right),$ (2)

$C_C={\left(u^{-{\theta }_2}+v^{-{\theta }_2}-1\right)}^{-1/{\theta }_2},$ (3)

$C_F=-\frac{1}{{\theta }_3}\text{log}\left(1+\frac{\left({\text{e}}^{-{\theta }_{3}u}-1\right)\left({\text{e}}^{-{\theta }_{3}v}-1\right)}{{\text{e}}^{-{\theta }_3}-1}\right).$ (4)

By adjusting the mixture weights, the mixture model can theoretically provide a more accurate approximation of the true dependence structure.

To quantify the dependence structure and tail synergy characteristics between NO₂ and PM_2.5, Kendall’s rank correlation coefficient $\tau$ , as well as upper tail dependence parameter ${\lambda }_U$ and lower tail dependence parameter ${\lambda }_L$ are employed [6] [26]:

$\tau =4{\displaystyle {\iint }_{I^2}C\left(u,v\right)\text{d}C\left(u,v\right)}-1,$ (5)

${\lambda }_U=2-\underset{t\to 1}{\text{lim}}\frac{1-C\left(t,t\right)}{1-t},$ (6)

${\lambda }_L=\underset{t\to 0}{\text{lim}}\frac{C\left(t,t\right)}{t}.$ (7)

Kendall’s $\tau$ reflects the degree of concordance in monotonic correlation between two variables [27]. A positive value indicates positive dependence, whereas a negative value indicates negative dependence. The tail dependence parameters quantify the probability of two variables significantly increasing or decreasing simultaneously within the same period. ${\lambda }_U$ measures the probability of joint extreme high values, while ${\lambda }_L$ measures the probability of joint extreme low values. Tail dependence parameters can reveal the synchronous fluctuation mechanisms of pollutants under heavy pollution or clean backgrounds, serving as important statistical indicators for characterizing extreme event risk coupling [28].

3.2. Marginal Distribution Estimation

Considering that pollutant concentration sequences typically exhibit heavy tails and non-normal characteristics, assuming parametric marginal distributions may introduce model misspecification bias. To enhance model robustness, this paper employs empirical distribution functions to estimate marginal distributions. For a sample size $n$ observation sequence $\left(X_1,Y_1\right),\cdots ,\left(X_n,Y_n\right)$ , the empirical distribution functions are defined as:

${\widehat{F}}_X\left(x\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^{n}I\left(X_i\le x\right),{\widehat{F}}_Y\left(y\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^{n}I\left(Y_i\le y\right),$(8)

where $I(\cdot )$ denotes the indicator function. Subsequently, rank transformation is used to convert the original observation sequence into pseudo-observations $\left(U_i,V_i\right)$ :

$U_i=\frac{\text{Rank}\left(X_i\right)}{n+1},V_i=\frac{\text{Rank}\left(Y_i\right)}{n+1}.$ (9)

3.3. Parameter Estimation

Maximum pseudo-likelihood estimation is widely used to estimate copula parameters [29]. This study employs maximum pseudo-likelihood estimation to solve the parameters of the mixture copula. The pseudo-log-likelihood function is defined as:

$\log L\left(\varphi \right)={\displaystyle \sum }_{i=1}^n\log \left({\alpha }_1{c}_G\left(u_i,v_i;{\theta }_1\right)+{\alpha }_2{c}_C\left(u_i,v_i;{\theta }_2\right)+{\alpha }_3{c}_F\left(u_i,v_i;{\theta }_3\right)\right),$ (10)

where $\varphi =\left({\alpha }_1,{\alpha }_2,{\theta }_1,{\theta }_2,{\theta }_3\right)$ is the parameter vector of the mixture copula, and $c_G,c_C,c_F$ are the probability density functions of Gumbel, Clayton and Frank copulas, respectively. Under parameter constraints, the interior point method is used to numerically maximize the objective function.

3.4. Model Evaluation and Goodness-of-Fit Testing

To systematically evaluate model performance, this paper constructs a multidimensional evaluation system, including the Akaike Information Criterion (AIC), root mean square error (RMSE), and Kolmogorov-Smirnov (KS) statistic [30] [31].

AIC balances fitting goodness and model complexity, defined as:

$\text{AIC}=2k-2\log L,$ (11)

where $k$ is the number of model parameters, and $\log L$ is the log-likelihood value.

Based on the pseudo-observations $\left(U_i,V_i\right)$ , the empirical copula is constructed as:

${\widehat{C}}_n\left(u,v\right)=\frac{1}{n}{\displaystyle \sum }_{i=1}^{n}I\left(U_i\le u,V_i\le v\right).$(12)

RMSE calculates the differences between the empirical copula and theoretical copula at all observation points:

$\text{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left({\widehat{C}}_n\left(u_i,v_i\right)-C\left(u_i,v_i\right)\right)}^2}.$(13)

The KS statistic tests the maximum deviation between empirical and model distributions:

$\text{KS}=\underset{i}{\text{max}}\left|{\widehat{C}}_n\left(u_i,v_i\right)-C\left(u_i,v_i\right)\right|.$(14)

To test the significance of the statistics, a parametric bootstrap procedure is implemented [32]. Simulated samples are generated based on estimated parameters, and RMSE and KS statistics are recalculated for each group of simulated samples, yielding $B$ simulated statistics. The p-value of the observed statistic in the simulated distribution is calculated using the +1 correction method to improve sample robustness:

$p=\frac{{{\displaystyle \sum }}_{b=1}^{B}I\left(T_b\ge T_{\text{obs}}\right)+1}{B+1},$ (15)

where $T_b$ is the statistic from the b-th simulation, and $T_{\text{obs}}$ is the observed statistic. A model is considered statistically adequate if the p-value exceeds the chosen significance level.

3.5. Clustering Analysis

Clustering analysis is a classical unsupervised learning method aimed at automatically grouping samples based on their similarity structure without prior knowledge [33]. As an important analytical tool in environmental science, clustering methods are widely applied to identify potential patterns and heterogeneity in complex observational datasets [34]. To distinguish differences among monitoring stations in terms of pollutant concentration characteristics and dependence structures, this paper adopts a clustering method based on statistical feature fingerprints. This method transforms raw observation sequences into low-dimensional feature vectors and achieves station grouping through traditional clustering algorithms. According to the type of pollution characteristics, clustering methods can generally be classified into shape-based, feature-based and model-based categories [35]. Given the obvious non-stationarity and extreme value sensitivity of pollutant data, this paper adopts a feature-based method, condensing high-dimensional observational data into interpretable statistical indicators for clustering analysis.

3.5.1. Statistical Feature Extraction

To enhance the robustness and interpretability of the analysis, the raw daily concentration series are not directly used. Instead, a statistical feature vector is constructed for each monitoring station, treating NO₂ and PM_2.5 as a coupled system for joint characterization. The feature vector consists of two categories:

1) Basic statistical features (12 dimensions). Mean, median, standard deviation, interquartile range (IQR), skewness and kurtosis are calculated separately for NO₂ and PM_2.5. These indicators comprehensively describe the average concentration levels, fluctuation degrees and distributional shape of pollutants. The IQR serves as a robust measure of variability that effectively resists extreme pollution values, while skewness and kurtosis reflect the asymmetry and tail thickness characteristics of the data distribution.

2) Interaction relationship features (6 dimensions). Pearson’s correlation coefficient between NO₂ and PM_2.5, linear regression slope, median concentration ratio, and high-concentration co-occurrence proportion (proportion where both exceed their respective 75th quantiles). These features quantify the linear dependence strength, relative growth relationships and extreme event synergy probabilities between pollutants, thereby facilitating the distinction of stations dominated by traffic sources from those dominated by regional transport.

3.5.2. Feature Standardization and Dimensionality Reduction

To eliminate differences in scale and mitigate potential collinearity among features, grouped Z-score standardization is applied to the original feature matrix. Subsequently, principal component analysis (PCA) is used to reduce the dimension of the standardized matrix, retaining principal components that collectively explain at least 90% of the cumulative variance [36]. The resulting low-dimensional score matrix obtained from PCA serves as the input for the clustering algorithm, reducing computational complexity while preserving the core structural information embedded in the original features.

3.5.3. Clustering Algorithm and Determination of the Optimal Number of Clusters

In the PCA-reduced space, the K-means algorithm is employed to cluster the monitoring stations [37]. To avoid local optima, 100 random initializations are set for repeated runs and the squared Euclidean distance is used as the similarity metric. Considering the potential spatial differentiation patterns within the study region, the number of clusters is pre-specified within the range $K=2$ to 5. The optimal number of clusters is determined based on the silhouette coefficient, which comprehensively evaluates within-cluster compactness and between-cluster separation. The $K$ value with the highest silhouette coefficient is selected as the optimal number of clusters.

4. Results and Discussion

4.1. Overall Joint Distribution Characteristics

To establish a foundational understanding of the joint distribution between NO₂ and PM_2.5 concentrations in Beijing, an overall analysis was conducted using daily average observations from 31 monitoring stations during 2023-2024. The overall sample results are shown in Table 1. During the study period, the average concentration of NO₂ is 25.704 μg/m³, whereas that of PM_2.5 is 33.456 μg/m³. The standard deviation of PM_2.5 (28.167 μg/m³) is substantially higher than that of NO₂ (13.695 μg/m³), indicating that PM_2.5 exhibits stronger volatility and instability. In terms of extreme values, the maximum value of PM_2.5 is 166.461 μg/m³, reflecting the occurrence of severe pollution episodes during the observation period. From a distributional perspective, the skewness of NO₂ and PM_2.5 is 0.957 and 1.816, respectively, and the corresponding kurtosis is 3.501 and 6.844. Both pollutants exhibit obvious right-skewed and heavy-tailed characteristics, with PM_2.5 displaying a markedly heavier tail. This suggests that extreme high-pollution events have a significant impact on the overall statistical structure.

Table 1. Descriptive statistics of NO₂ and PM_2.5 (μg/m³).

Variable	Mean	Min	Max	Std. Dev.	Skewness	Kurtosis
NO2	25.704	4.853	72.250	13.695	0.957	3.501
PM2.5	33.456	1.658	166.461	28.167	1.816	6.844

Figure 2 shows the scatter plot between the two variables, showing an overall positive correlation between them. Data points are densely clustered in the low-concentration range, exhibiting an approximately linear trend. Observations in the high-concentration range become sparse and more dispersed, without an evident clustering pattern.

Figure 2. Scatter plot of NO₂ and PM_2.5.

To further quantify the nonlinear and asymmetric joint distribution characteristics between NO₂ and PM_2.5 concentration sequences, based on the above joint samples, we transformed the marginal distributions using the empirical distribution function obtained from Equation (8), and then fitted Gumbel, Clayton, Frank copulas and their mixture copula model, respectively.

Model selection was conducted using the AIC, RMSE and the KS test. The model evaluation results are shown in Table 2. The mixture copula model has an AIC value of −382.844 and RMSE of 0.007, both optimal among the four. Moreover, it passes the KS test, with a KS statistic of 0.020 and a p-value of 0.276, indicating no statistically significant difference between the fitted and empirical distributions at the 0.05 significance level. Statistically, the mixture copula demonstrates the best overall goodness-of-fit.

Table 2. Comparison of copula model fitting goodness for overall data.

Copula	AIC	RMSE	$\tau$	${\lambda }_U$	${\lambda }_L$
Mixture	−382.844	0.007	0.451	0.265	0.085
Gumbel	−370.890	0.009	0.422	0.507	0
Clayton	−295.595	0.023	0.346	0	0.520
Frank	−373.525	0.011	0.451	0	0

To visually verify the ability of each model to characterize the data dependence structure, the empirical joint probability density function is compared with the theoretical densities implied by the four copula models (Figure 3 and Figure 4). From the figures, it is clear that a single copula model has obvious limitations in capturing dependence structures. Gumbel copula tends to overemphasize upper tail dependence, Clayton copula is more sensitive to lower tail dependence, and Frank copula captures only symmetric dependence. In contrast, the theoretical density of the mixture copula has the highest correlation with the empirical density, not only accurately depicting the dependency patterns within the main data range, but also effectively capturing the asymmetric features at the tails.

Figure 3. Empirical joint probability density.

Figure 4. Theoretical joint probability density: (a) Mixture; (b) Gumbel; (c) Clayton; (d) Frank.

Figure 5. Q-Q plot of the mixture copula.

After determining the mixture copula as the optimal model, its quantile fitting performance was further examined using a Q-Q plot comparing theoretical and empirical quantiles (Figure 5). From Figure 5, it can be seen that sample points are basically distributed around the diagonal, with no obvious deviation even in extreme quantile intervals close to 0 and 1. It indicates that the mixture copula model not only fits well overall but also accurately characterizes the dependence structure of NO₂ and PM_2.5 under extreme concentrations.

Substituting the NO₂ and PM_2.5 data into Equations (1)-(4) and Equation (10) yielded the parameters for the four copula models, with parameter estimation results shown in Table 3. For the mixture copula, the weight parameters are ${\alpha }_1=0.523$ , ${\alpha }_2=0.131$ , ${\alpha }_3=0.346$ . These results indicate that the model is upper-tail-dominant dependence structure of Gumbel copula, while incorporating lower-tail-dominant dependence structure from Clayton copula and symmetric dependence from Frank copula. This weighted combination produces an asymmetric dependence structure more consistent with reality. The Kendall’s $\tau$ of the mixture copula is 0.451, with an upper tail dependence parameter ${\lambda }_U$ of 0.265 and a lower tail dependence parameter ${\lambda }_L$ of 0.085. This indicates that a moderate positive association between NO₂ and PM_2.5, with a substantially higher probability of joint extreme high concentrations than joint extreme low concentrations. The dependence structure can therefore be characterized as being dominated by high-pollution synergy.

Table 3. Parameter estimates of the four copula models.

Copula	${\theta }_1$	${\theta }_2$	${\theta }_3$	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$
Mixture	1.730	1.625	5.564	0.523	0.131	0.346
Gumbel	1.730	–	–	–	–	–
Clayton	–	1.059	–	–	–	–
Frank	–	–	4.914	–	–	–

In summary, the mixture copula not only comprehensively captures the multiple dependence patterns between NO₂ and PM_2.5 concentrations but also performs best across all fitting goodness indicators, making it suitable for describing complex dependence structures in such environmental pollutant data.

4.2. Spatial Clustering Analysis of Stations

To investigate the influence of spatial heterogeneity on the dependence structure between NO₂ and PM_2.5, the K-means algorithm was applied to cluster the 31 air quality monitoring stations. The results indicate that the silhouette coefficient reaches its maximum value (0.692) when $K=2$ , suggesting that partitioning the stations into two groups yields the optimal clustering performance. This configuration effectively distinguishes the core differences in pollutant concentration characteristics and dependency relationships across regions. The final clustering results are to divide the 31 stations into two groups, with spatial distribution as shown in Figure 6, highly consistent with Beijing’s “urban-suburban” geographical pattern. In subsequent analysis, the first group from the clustering results is uniformly referred to as the “urban-biased group,” and the second group as the “suburban-biased group.”

Figure 6. Spatial distribution of clustering.

The urban-biased group consists of 18 stations, including Dongcheng Dongsi, Dongcheng Tiantan, Fengtai Yungang, Fengtai Xiaotun, Yizhuang Development Zone, Jingdongnan Regional Point, Daxing Jiugong, Daxing Huangcun, Fangshan Liangxiang, Chaoyang Nongzhanguan, Chaoyang Aoti Center, Haidian Wanliu, Shijingshan Gucheng, Shijingshan Laoshan, Xicheng Wanshouxigong, Xicheng Guanyuan, Tongzhou Dongguan, Tongzhou Yongshun. These stations are mainly located within central urban districts and southeastern near-suburban areas. Their locations are close to major urban traffic arteries, commercial areas and densely populated regions, exhibiting typical urban pollution characteristics. The concentration levels of NO₂ and PM_2.5 are generally high, with average values of 30.63 μg/m³ and 35.88 μg/m³, respectively. The Pearson’s correlation coefficient is 0.618, showing strong correlation. This pattern reflects the synergistic pollution characteristics dominated by local emission sources, particularly traffic-related emissions in densely populated areas with high human activity.

The suburban-biased group consists of 13 stations, including Dingling Control Point, Miyun New City, Miyun Town, Pinggu New City, Pinggu Town, Yanqing Xiandu, Yanqing Shiheying, Huairou Town, Fangshan Yanshan, Changping Town, Haidian Sijiqing, Mentougou Shuangyu, and Shunyi Beixiaoying. These stations are mainly distributed in the northwestern and northeastern distant suburbs and ecological conservation zones, far from the urban core, with relatively less influence from local emissions. The overall pollution level of this group of stations is relatively low, with the average values of NO₂ and PM_2.5 being 18.89 μg/m³ and 30.10 μg/m³, respectively. Notably, NO₂ levels are substantially lower than those in the urban-biased group, indicating reduced anthropogenic emission intensity. The Pearson’s correlation coefficient is 0.568, slightly lower than that of the urban-biased group, implying that PM_2.5 sources in these areas may be more complex and relatively more influenced by regional transport and secondary generation processes.

The above clustering results were basically consistent with Beijing’s actual spatial pattern, indicating that the clustering method effectively captured emission structure differences. The key difference between the two groups lies in the obvious greater decline in NO₂ than PM_2.5, indicating that atmospheric pollutant concentrations in Beijing exhibited obvious spatial heterogeneity. The urban-biased group stations are dominated by local primary emission sources (traffic, industry, domestic, etc.), while suburban-biased group stations more reflect regional transport and natural background characteristics, with a higher proportion of regional transport and secondary generation contributions to PM_2.5 in suburbs. Meanwhile, the correlation differences between the two groups are limited, indicating that spatial location primarily affects pollutant concentration levels, with limited impact on the overall form of the dependence structure.

4.3. Seasonal Analysis of Pollutant Concentrations

This subsection further takes a temporal dimension into consideration, analyzing the impact of seasonal changes on the dependence structure between pollutants to further elucidate the dynamic association characteristics of pollutants.

Figure 7 shows the seasonal variations in monthly average concentrations of NO₂ and PM_2.5 in Beijing from 2023-2024. A highly consistent seasonal cycle is observed for both pollutants: concentrations decrease synchronously during spring and summer, reaching their lowest levels in summer. In contrast, concentrations increase markedly in autumn and winter, peaking during winter months. According to the official heating season classification in Beijing, the heating season is from November 15 to March 15 of the following year, and the non-heating season is from March 16 to November 14. The overall concentration levels of both pollutants during the heating season are significantly higher than those during the non-heating season. This phenomenon reflects that anthropogenic emission intensity such as coal burning significantly increased during the heating season, and meteorological factors such as stable winter atmospheric boundary layers and poor pollutant diffusion conditions jointly lead to enhanced pollution accumulation effects. The seasonal differences in pollutant concentrations suggest that the dependence structure between NO₂ and PM_2.5 might undergo structural reconfiguration under varying emission intensities and meteorological backgrounds. Therefore, to deeply analyze their joint dependence structure, it is necessary to break through the limitations of full-year mixed analysis and stratify the data by heating and non-heating seasons to capture their potential seasonal pattern transition characteristics.

Figure 7. Seasonal variations in monthly average concentrations of NO₂ and PM_2.5 in Beijing from 2023-2024: (a) 2023; (b) 2024.

4.4. Dependence Structure Modeling from a Spatiotemporal Heterogeneity Perspective

Based on the preceding clustering and seasonal analyses, there are obvious regional and seasonal differences between NO₂ and PM_2.5 concentrations. Accordingly, a two-dimensional analytical framework integrating spatial grouping and seasonal stratification was constructed, yielding four distinct scenarios: urban-biased-heating season, suburban-biased-heating season, urban-biased-non-heating season, and suburban-biased-non-heating season. In each scenario, Gumbel, Clayton, Frank and mixture copula models were fitted separately, enabling a systematic comparison of model performance and dependence characteristics under different spatiotemporal conditions.

Table 4 summarizes the optimal copula model and parameter estimation results for each scenario. At the significance level $p=0.05$ , all scenario optimal models pass the KS test, indicating no significant difference between the fitted and empirical distributions, confirming the reliability of the models.

Table 4. Optimal fitting models for the dependence structure between NO₂ and PM_2.5 under different spatiotemporal scenarios.

Region	Period	Copula	$\tau$	${\lambda }_U$	${\lambda }_L$
Urban-biased group	Heating season	Mixture ($0.097C_G+0.506C_C+0.397C_F$ )	0.558	0.063	0.380
Urban-biased group	Non-heating season	Gumbel	0.367	0.449	0
Suburban-biased group	Heating season	Mixture ($0.602C_C+0.398C_F$ )	0.557	0	0.450
Suburban-biased group	Non-heating season	Gumbel	0.384	0.467	0

During the heating season, the dependence structure between NO₂ and PM_2.5 in each region exhibits high complexity and asymmetry. The optimal models are all mixture copulas, indicating that a single copula model is insufficient to adequately capture the composite joint distribution patterns. Kendall’s $\tau$ exceeds 0.5 in both cases, reflecting moderate positive dependence between pollutants. Tail dependence analysis further revealed that lower tail dependence (${\lambda }_L=0.3\text{-}0.5$ ) is stronger than upper tail dependence (${\lambda }_L<0.1$ ). This characteristic is jointly driven by the emission patterns during the heating season, boundary-layer dynamics, and secondary formation [38]-[40]. Specifically, coal-fired heating induces a surge in PM_2.5 emissions, accompanied by simultaneous increases in NO₂ precursors, resulting in highly coupled emission sources [38]. This elevates the overall pollution level and shifts high-concentration events from extreme occurrences to relatively common conditions, thereby weakening the statistical significance of upper tail dependence. Moreover, the heating season is characterized by a shallow boundary layer and frequent stagnant conditions, which hinder pollutant diffusion [39]. Both pollutants can decline simultaneously only under favorable dispersion conditions, such as strong winds or precipitation. Therefore, synchronized fluctuations are more likely to occur in the lower quantile range, strengthening the lower tail dependence. In addition, low-temperature and stable environments facilitate the conversion of NO₂ into nitrate aerosols, establishing a chemical source-sink relationship. In this process, NO₂ serves as the precursor and PM_2.5 as the nitrate reservoir [40], which enhances the overall positive dependence. Consequently, the heating season dependence structure can be characterized by moderate positive dependence with a lower-tail-dominant pattern.

In contrast, the dependence structure during the non-heating season is relatively simple. For both spatial groups, Gumbel copula provides the optimal fit. This indicates that a single copula model is sufficient for modeling and interpretation, whereas the additional parameters of the mixture copula enhance flexibility, also significantly increase the estimation complexity, thereby resulting in a decrease in fitting efficiency. Kendall’s $\tau$ ranges between 0.3 and 0.4, indicating moderate positive dependence. The upper tail dependence parameters exceed 0.4, while lower tail dependence parameters are all 0. During the non-heating season, the atmospheric boundary-layer height increases and convective activity intensifies, facilitating pollutant dispersion. Synergistic accumulation at low concentrations is unlikely to occur, and lower tail dependence is weak [10]. Conversely, stagnant meteorological conditions under high temperature and humidity markedly enhance atmospheric oxidizability, promoting photochemical smog formation and secondary aerosol production [41] [42]. This drives the simultaneous increase of NO₂ and PM_2.5, resulting in a strong upper tail dependence.

Within the same season, the dependence patterns of urban-biased and suburban-biased groups stations are similar, but subtle differences exist. For instance, during the heating season, the lower tail dependence in the suburban-biased group is slightly stronger than in the urban-biased group. This difference might reflect the relatively weaker influence of direct anthropogenic emissions in suburban environments, making them more sensitive to synchronized pollutant reductions driven by natural deposition or regional transport dilution processes [10]. Overall, seasonal factors exerted a markedly stronger regulatory influence on the dependence structure between NO₂ and PM_2.5 than spatial differences. Spatial location primarily affects absolute concentration levels, whereas its structural impact on the coupling mechanisms remained limited.

The hierarchical analysis demonstrates that the dependence structure between NO₂ and PM_2.5 in Beijing is not static but exhibits significant spatiotemporal heterogeneity. The overall tail dominance pattern is decomposed into lower tail dominance during the heating season and upper tail dominance during the non-heating season after stratification. This structural transformation arises from multiple mechanisms, including data aggregation effects (overall averaging obscures local patterns), differentiation of emission sources (urban local emissions versus suburban regional transport), and seasonal meteorological regulation (winter accumulation versus summer photochemical generation). These findings indicate the critical role of a spatiotemporal analytical framework in revealing compound pollution dynamics.

4.5. Limitations

This study has several limitations. First, by treating the heating season and non-heating season as holistic analytical units, the model does not fully account for non-stationarity within seasons. This simplification may obscure dynamic variations at finer temporal scales. Second, meteorological factors (such as wind speed and temperature) were not incorporated as covariates, potentially limiting a deeper mechanistic interpretation of the dependence structure. Finally, the estimated dependence structures reflect statistical associations rather than causal relationships. The joint fluctuations of NO₂ and PM_2.5 result from multiple interacting factors, including emissions, meteorology, and chemical transformations, and copula parameters cannot be interpreted as evidence of causal mechanisms.

5. Conclusions

This paper employs an analytical framework integrating overall modeling, spatial clustering and seasonal stratification. Based on daily average NO₂ and PM_2.5 data from 31 air quality monitoring stations in Beijing during 2023-2024, it systematically reveals the spatiotemporal heterogeneity of pollutant joint dependence structures. The main conclusions are as follows:

First, at the overall scale, NO₂ and PM_2.5 in Beijing exhibit a moderate positive dependence characterized predominantly by asymmetric upper-tail-dominant dependence structure. This indicates a joint amplification effect in high-pollution scenarios. Second, at the spatial scale, stations are clustered into urban-biased and suburban-biased groups, with obvious gradient differences in concentration levels (especially NO₂ spatial stratification), but the dependence structure remains relatively stable. Spatial location primarily regulates pollution intensity rather than coupling mechanisms between pollutants. Third, at the temporal scale, both NO₂ and PM_2.5 show clear seasonal cycles. Concentration levels of two pollutants increase significantly during the heating season and decline markedly during the non-heating season, emphasizing the critical roles of meteorological conditions and anthropogenic emissions in pollution accumulation. Fourth, the spatiotemporal heterogeneity analysis indicates that seasonal changes significantly restructure the dependence pattern. The heating season is characterized by strong overall positive dependence and lower-tail-dominant dependence, whereas the non-heating season exhibits upper-tail-dominant dependence. Seasonal effect appears stronger than spatial differences, indicating that the dependence relationship between pollutants has spatiotemporal heterogeneity.

From a policy perspective, local fitting yields substantially greater practical relevance than overall fitting. While the overall model provides a citywide baseline reference, it is difficult to support differentiated strategies. Spatiotemporal stratification directly points to targeted measures: urban-biased stations need to prioritize regulation of traffic and industrial sources to mitigate the risk of joint high-concentration events; suburban-biased stations require strengthening regional joint prevention and reduce the impact of pollution transmission; during the heating season, emission limits should be implemented, and during non-heating seasons, particular attention should be paid to controlling NO₂ as a key ozone precursor to suppress photochemical synergistic amplification. Such a refined strategy is far more effective than uniform, citywide interventions.

Overall, this study demonstrates that the relationship between NO₂ and PM_2.5 in Beijing is not a stable linear correlation, but undergoes structural reconfiguration under varying emission and meteorological conditions. These findings provide important insights for understanding synergistic mechanisms of compound pollution, assessing joint risks of extreme pollution events, and formulating seasonally differentiated control strategies. Furthermore, the proposed framework establishes a methodological foundation for extending multivariate copula models in future studies.

Acknowledgements

Sincere thanks to the members of JAMP for their professional performance, and special thanks to managing editor Hellen XU for a rare attitude of high quality.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	He, C., Pan, F. and Yan, Y. (2011) Is Economic Transition Harmful to China’s Urban Environment? Evidence from Industrial Air Pollution in Chinese Cities. Urban Studies, 49, 1767-1790.[CrossRef]
[2]	Ebenstein, A., Fan, M., Greenstone, M., He, G., Yin, P. and Zhou, M. (2015) Growth, Pollution, and Life Expectancy: China from 1991-2012. American Economic Review, 105, 226-231.[CrossRef]
[3]	Li, G., Li, L., Liu, D., Qin, J. and Zhu, H. (2021) Effect of PM2.5 Pollution on Perinatal Mortality in China. Scientific Reports, 11, Article No. 7596.[CrossRef] [PubMed]
[4]	Li, S., Zhou, C. and Qu, M. (2023) Spatiotemporal Variations and Mechanism of PM2.5 Pollution in Urban Area: The Case of Guiyang, Guizhou, China. Journal of Environmental Management, 341, Article ID: 118030.[CrossRef] [PubMed]
[5]	Lim, S.S., Vos, T., Flaxman, A.D., Danaei, G., Shibuya, K., Adair-Rohani, H., et al. (2012) A Comparative Risk Assessment of Burden of Disease and Injury Attributable to 67 Risk Factors and Risk Factor Clusters in 21 Regions, 1990-2010. The Lancet, 380, 2224-2260.
[6]	Nelsen, R.B. (2006) An Introduction to Copulas. Springer.
[7]	Masseran, N. (2021) Modeling the Characteristics of Unhealthy Air Pollution Events: A Copula Approach. International Journal of Environmental Research and Public Health, 18, Article No. 8751.[CrossRef] [PubMed]
[8]	Kumar, P. (2018) Copula Functions and Applications in Engineering. In: Deep, K., Jain, M. and Salhi, S., Eds., Logistics, Supply Chain and Financial Predictive Analytics: Theory and Practices, Springer, 195-209.[CrossRef]
[9]	Oppenheimer, M., Little, C.M. and Cooke, R.M. (2016) Expert Judgement and Uncertainty Quantification for Climate Change. Nature Climate Change, 6, 445-451.[CrossRef]
[10]	Abbasi, M.T., Alesheikh, A.A., Jafari, A. and Lotfata, A. (2024) Spatial and Temporal Patterns of Urban Air Pollution in Tehran with a Focus on PM2.5 and Associated Pollutants. Scientific Reports, 14, Article No. 25150.[CrossRef] [PubMed]
[11]	Zhang, J., Li, Y., Liu, C., Wu, B. and Shi, K. (2022) A Study of Cross-Correlations between PM2.5 and O3 Based on Copula and Multifractal Methods. Physica A: Statistical Mechanics and Its Applications, 589, Article ID: 126651.[CrossRef]
[12]	Sohrabian, B. (2021) Geostatistical Prediction through Convex Combination of Archimedean Copulas. Spatial Statistics, 41, Article ID: 100488.[CrossRef]
[13]	Masseran, N. and Hussain, S.I. (2020) Copula Modelling on the Dynamic Dependence Structure of Multiple Air Pollutant Variables. Mathematics, 8, Article No. 1910.[CrossRef]
[14]	Sak, H., Yang, G., Li, B. and Li, W. (2017) A Copula-Based Model for Air Pollution Portfolio Risk and Its Efficient Simulation. Stochastic Environmental Research and Risk Assessment, 31, 2607-2616.[CrossRef]
[15]	Chan, R.K.S. and So, M.K.P. (2018) Multivariate Modelling of Spatial Extremes Based on Copulas. Journal of Statistical Computation and Simulation, 88, 2404-2424.[CrossRef]
[16]	Qian, L., Wang, X. and Wang, Z. (2020) Modeling the Dependence Pattern between Two Precipitation Variables Using a Coupled Copula. Environmental Earth Sciences, 79, Article No. 486.[CrossRef]
[17]	Veerman, R. (2022) Modeling Dependence between Financial Markets: A Copula Approach, Master’s Thesis, University of Amsterdam.
[18]	Gontara, E. and Chebana, F. (2025) Mixture Copula Parameter Estimation with Metaheuristic Algorithms, Comparative Study under Hydrological Context. Stochastic Environmental Research and Risk Assessment, 39, 1307-1326.[CrossRef]
[19]	Arakelian, V. and Karlis, D. (2014) Clustering Dependencies via Mixtures of Copulas. Communications in Statistics—Simulation and Computation, 43, 1644-1661.[CrossRef]
[20]	Chi, Y., Ren, Y., Xu, C. and Zhan, Y. (2024) The Spatial Distribution Mechanism of PM2.5 and NO2 on the Eastern Coast of China. Environmental Pollution, 342, Article ID: 123122.[CrossRef] [PubMed]
[21]	Gómez-Carracedo, M.P., Andrade, J.M., López-Mahía, P., Muniategui, S. and Prada, D. (2014) A Practical Comparison of Single and Multiple Imputation Methods to Handle Complex Missing Data in Air Quality Datasets. Chemometrics and Intelligent Laboratory Systems, 134, 23-33.[CrossRef]
[22]	Umar, N. and Gray, A. (2023) Comparing Single and Multiple Imputation Approaches for Missing Values in Univariate and Multivariate Water Level Data. Water, 15, Article No. 1519.[CrossRef]
[23]	Junninen, H., Niska, H., Tuppurainen, K., Ruuskanen, J. and Kolehmainen, M. (2004) Methods for Imputation of Missing Values in Air Quality Data Sets. Atmospheric Environment, 38, 2895-2907.[CrossRef]
[24]	Durante, F. and Sempi, C. (2015) Principles of Copula Theory. CRC Press.[CrossRef]
[25]	Flores, M.U., de Amo, E., Durante, A.F. and Sánchez, J.F. (2017) Copulas and Dependence Models with Applications. Springer.
[26]	Talbi, M., de Peretti, C. and Belkacem, L. (2020) Dynamics and Causality in Distribution between Spot and Future Precious Metals: A Copula Approach. Resources Policy, 66, Article ID: 101645.[CrossRef]
[27]	Han, S., Qiao, Y., Yan, J., Liu, Y., Li, L. and Wang, Z. (2019) Mid-to-Long Term Wind and Photovoltaic Power Generation Prediction Based on Copula Function and Long Short Term Memory Network. Applied Energy, 239, 181-191.[CrossRef]
[28]	Zhang, Y., Kim, C.W., Beer, M., Dai, H. and Soares, C.G. (2018) Modeling Multivariate Ocean Data Using Asymmetric Copulas. Coastal Engineering, 135, 91-111.[CrossRef]
[29]	Qu, L. and Lu, Y. (2019) Copula Density Estimation by Finite Mixture of Parametric Copula Densities. Communications in Statistics - Simulation and Computation, 50, 3315-3337.[CrossRef]
[30]	Zucchini, W. (2000) An Introduction to Model Selection. Journal of Mathematical Psychology, 44, 41-61.[CrossRef] [PubMed]
[31]	Joe, H. (2014) Dependence Modeling with Copulas. CRC Press.
[32]	Genest, C., Rémillard, B. and Beaudoin, D. (2009) Goodness-of-Fit Tests for Copulas: A Review and a Power Study. Insurance: Mathematics and Economics, 44, 199-213.[CrossRef]
[33]	Jain, A.K., Murty, M.N. and Flynn, P.J. (1999) Data Clustering: A Review. ACM Computing Surveys, 31, 264-323.[CrossRef]
[34]	D’Urso, P., De Giovanni, L. and Massari, R. (2015) Time Series Clustering by a Robust Autoregressive Metric with Application to Air Pollution. Chemometrics and Intelligent Laboratory Systems, 141, 107-124.[CrossRef]
[35]	Liao, T.W. (2005) Clustering of Time Series Data—A Survey. Pattern Recognition, 38, 1857-1874.[CrossRef]
[36]	Lu, W.Z., He, H.D. and Dong, L.Y. (2011) Performance Assessment of Air Quality Monitoring Networks Using Principal Component Analysis and Cluster Analysis. Building and Environment, 46, 577-583.[CrossRef]
[37]	Ahmed, M., Seraj, R. and Islam, S.M.S. (2020) The K-Means Algorithm: A Comprehensive Survey and Performance Evaluation. Electronics, 9, Article No. 1295.[CrossRef]
[38]	Wan, Z., Yang, C., Wang, X., Xue, Y., Zhao, J., Cui, J., et al. (2024) Spatial and Temporal Differentiation of Air Quality and Its Influence Factors in 16 Cities in Shandong Province from 2019 to 2020. RSC Sustainability, 2, 1528-1542.[CrossRef]
[39]	Zhong, J., Zhang, X., Dong, Y., Wang, Y., Liu, C., Wang, J., et al. (2018) Feedback Effects of Boundary-Layer Meteorological Factors on Cumulative Explosive Growth of PM2.5 during Winter Heavy Pollution Episodes in Beijing from 2013 to 2016. Atmospheric Chemistry and Physics, 18, 247-258.[CrossRef]
[40]	Huang, R.J., Zhang, Y., Bozzetti, C., Ho, K.F., Cao, J.J., Han, Y., et al. (2014) High Secondary Aerosol Contribution to Particulate Pollution during Haze Events in China. Nature, 514, 218-222.[CrossRef] [PubMed]
[41]	Zhao, D., Liu, G., Xin, J., Quan, J., Wang, Y., Wang, X., et al. (2020) Haze Pollution under a High Atmospheric Oxidization Capacity in Summer in Beijing: Insights into Formation Mechanism of Atmospheric Physicochemical Processes. Atmospheric Chemistry and Physics, 20, 4575-4592.[CrossRef]
[42]	Wang, D., Zhou, B., Fu, Q., Zhao, Q., Zhang, Q., Chen, J., et al. (2016) Intense Secondary Aerosol Formation Due to Strong Atmospheric Photochemical Reactions in Summer: Observations at a Rural Site in Eastern Yangtze River Delta of China. Science of the Total Environment, 571, 1454-1466.[CrossRef] [PubMed]