Robust Machine Learning Algorithmic Rules for Detecting Air Pollution in the Lower Parts of the Atmosphere

1.2.1 Complexity of SDG Interactions

For an intuition into the mutual impact of human–nature interactions, consider the inherent correlations among different aspects of the 17 SDGs (United Nations, 2015). One aspect that relates to each SDG is “environmental pollution”. Thus, conducting air quality studies and spatial analyses can play an important role in planning and laying down policies for land use, infrastructural and business development, which hinge on different SDGs. For instance, if pollution hotspots in southern China (a major trading port and manufacturing hub) can be identified, preventive measures and stricter pollution abatement policies can be implemented in concerned cities or counties, so that overall pollution levels can be effectively reduced, ensuring healthy lives and promoting well–being for all citizens and visitors (SDG #3).

Air quality and climatic conditions mutually impact each other: climate change can affect air quality, and certain air pollutants can affect climate change (Feng et al., 2019) (SDG #13). Other SDGs that are closely associated with air quality include SDG #11 (Sustainable Cities and Communities), the attainment of which relates to reducing air pollution and making cities and human settlements inclusive, safe, resilient and sustainable. Smart and sustainable city development shows care and a focus on maintaining good air quality within neighbourhood levels and within various spatial scales. Focusing on a specific region, such as southern China, contributes to SDG #17 (Partnerships), which facilitates interdisciplinarity across geographical regions, legislations, infrastructural setups, and promotes citizen engagement.

1.2.2 Addressing the Complexity of SDG interactions

Recognising the complexity stated in Section 1.2.1, we adopt the mathematical concept of “non-orthogonality” of SDGs (Mwitondi, Munyakazi and Gatsheni, 2020). Different methods have been developed to track, monitor and model environmental pollution, including automated data-driven tools to capture, model and track environmental variations (Mwitondi, Munyakazi and Gatsheni, 2020). The literature on environmental pollution is awash with applications of sophisticated statistical and machine learning methods (Liu et al., 2022a; Pan, Harrou and Sun, 2023). Data–intensive, machine learning approaches have been used in the remote processing of ground-level concentrations of atmospheric pollutants such as ${\text{NO}}_2$, PM_2.5 and ozones (Chi et al., 2022; Du et al., 2022; Xu et al., 2018). Clustering methods have been applied in quantifying pollution levels as well as in remote sensing (Hsu et al., 2023; Zhang and Yang, 2022). Other studies have applied temperature inversion methods (Feng, Wei and Wang, 2020), partial differential equations with appropriate initial and boundary conditions (Shafiev, 2024), as well as the use of remotely sensed datasets (Lin et al., 2020; Mak et al., 2018). A study focusing on surface and sea surface pressure, geopotential height, temperature, relative humidity, wind field, and vertical velocity, established an association between regional weather and climate events on one side and large-scale circulation anomalies on the other (Cai et al., 2020). Global studies on the relationship between economic development and pollution are also well documented (Yan et al., 2024). The main challenge of the foregoing studies revolves around model optimisation (Liu et al., 2022a; Vardoulakis et al., 2007; Zhang et al., 2022a), i.e., the need to develop “robust” methods to perform relevant tasks that can be replicated in a spatio-temporal context. Further, socio-demographic factors, coupled with economic, cultural, political, legislative and technical variations, directly affect pollution levels as well as mitigation efforts. Focusing on southern China was therefore particularly appealing due to the country’s role on both sides of the pollution spectrum, i.e., as a polluter and as a pollution mitigator (Liu et al., 2022b; Sun, 2016; Zheng and Kahn, 2017). Furthermore, most of the aforementioned studies have been carried out in China, hence associating results from a “robust” model with “acceptable metrics” such as SDG indicators and sharing the relevant datasets publicly would be a major step towards our understanding of SDG attainment (Mak and Lam, 2021; Mwitondi, Munyakazi and Gatsheni, 2020).

Our paper complements previous studies by adopting an interdisciplinary approach to spatio-temporal air pollution modelling, via the two algorithms in Section 2.2.4 that are embedded with the flexibility to run different techniques including, in this specific application, Principal Component Analysis (PCA), Cluster Analysis (Chapmann, 2017; Kogan, 2007), Correspondence Analysis (CA) (Hirschfeld, 1935), K-Means (Lloyd, 1957; MacQueen, 1967) and the Expectation-Maximisation (EM) algorithm (Dempster et al., 1977). Our main idea hinges on an interdisciplinary approach to both problem identification and solving, i.e., balancing the power of data, machine learning techniques, and underlying domain knowledge. We focus on the identification of pollution levels of different pollutants in the sampled area, then propose two algorithms for identifying spatio-temporal associations among key air pollutants. The novelty of this study derives from the role of the two algorithms, which hinges on their “robustness” in addressing data randomness and on highlighting the path towards addressing the “non-orthogonality” of SDGs, i.e., the factors that affect the extent and impact of air pollution within any geographical area are not confined to the geographical boundaries of that area. This can be better understood in the context of how our atmosphere is structured, because pollutants along the air can travel from one country to another within a short period of time, as well as across various vertical layers within the atmosphere. Table 1 describes the structured layers of our atmosphere and their relevance to humanity.

In this study, we focus on the tropospheric layer – the bottom layer of the earth’s atmosphere. It constitutes of about 75–80% of the atmospheric mass, with its temperature variation affected by height and time of day and/or year. Most of our terrestrial weather—clouds, rain, snow—derive from it, making it a particularly interesting research scope.

2.2.1 Computational Maximisation of the Log-Likelihood

Some aspects of data visualisation are presented in Section 3.1 to convert a complex scenario into an easy one, which is useful for readers with limited knowledge of statistics. The popularity of data visualisation across data science applications is well documented (Liang et al., 2023; Pika et al., 2021). However, data visualisation has its limitations, particularly when the interest lies on representing characteristics of high-dimensional data. Under such circumstances, density estimation methods, maximum likelihood and posterior estimations are the most appropriate methods.

It is reasonable to assume that the collated data variables are independent and identically distributed, i.e., they are mutually independent, with the same probability distribution. Let us denote the five pollutants by $K$ simple distributions

where ${\varphi }_k$ represent basic distributions and ${\pi }_k$ are group proportions or mixture weights. This implies that we can describe multimodal distributions of any dimension, which we would otherwise only clearly visualise in a 1, 2 or 3-dimensional space. Consider the components of a Gaussian Mixture Model, each with parameters ${\mu }_k$ and ${\mathrm{\Sigma }}_k$.

where $\theta =\{{\pi }_k,{\mu }_k,{\mathrm{\Sigma }}_k;k=1,2,\dots ,K–1,K\}$ are the free parameters of the distribution. Let us denote the pollution data in Table 2 by $\mathbf{\text{X}}=\{x_1,x_2,\dots x_N\}$ where $x_n,n=1,2,\dots ,N–1,N$ are independent and identically distributed from some unknown distribution $\varphi (x).$ It can be shown that Equation 2 leads to a set of dependent simultaneous equations that can be solved iteratively. To estimate the aggregate distribution representing density distributions such as those exhibited in Figure 5, we can initialise $K$ mixture components as follows:

Given ${\overline{x}}_k\text{and}{\sigma }_k,$ the estimated aggregating density in Equation 2 will depend on the weights $({\pi }_k)$ assigned to each of the $K$ components, as ${\pi }_k$ defines prior memberships to each–i.e., prior probability of the type of pollution. The maximum likelihood estimate of the free parameters ${\theta }_{\text{ML}}$ is given as follows:

where each of the likelihood terms is assumed to be a Gaussian Mixture Model density, with the log likelihood

For a single Gaussian model, the sum over $K$ vanishes, and Equation 5 reduces to

The parameters that maximise the log likelihood in Equation 5 are obtained by maximising the derivatives of the mean ($\mu$), variation ($\mathrm{\Sigma }$) and group proportions ($\pi$) with respect to $\theta$, as follows

Note that the expression in Equation 7, i.e., the necessary conditions for maximising Equation 5, are of the form

where $\theta =\{{\mu }_k,{\mathrm{\Sigma }}_k,{\pi }_k;k=1,2,\dots ,K\}$ are model parameters, $p(x_n|\theta )$ is the probability of data $x_n$, given $\theta ,$ hence

The expressions in the denominator in Equation 9 imply that the probability of data given parameters is proportional to the sum of all the $K$ components. It can be further shown that

The left-hand side of Equation 10 is a normalised probability vector that represents the probability that the data $\left(x_n\right)$ was generated by the $k^{\text{th}}$ mixture component, which is directly proportional to the likelihood on the right-hand side. Direct maximisation of the likelihood function of data in the form of a random sample is challenging, so the likelihood and parameters in Equations 1 to 10 can only be obtained by iteratively searching for parameters that maximise $\mathcal{L}.$

2.2.2 Maximum-Likelihood Parameter Estimation and Maximisation Algorithm

If we envision each observation as being characterised by a parametric finite mixture density, we can treat group membership as missing data and use an adapted version of the EM algorithm to estimate and maximise the parameters ${\mu }_k$ and ${\mathrm{\Sigma }}_k.$ Let $c_i=1,2,\mathrm{\dots },K$ be initial classes with an unobservable indicator variable

We can then compute class membership as the probability of class given data, as follows

If class membership, central tendency, and variation parameters were observable, they could be estimated as

where ${\widehat{\pi }}_k$ is the estimated proportion of data points in class $k,{\widehat{\mu }}_k$ is the mean within that class and ${\widehat{\sigma }}_k^2$ is the variation in that category. However, the parameters are not observable, but we can estimate them from data, at each E step

The estimated parameters are then maximised at the $M$ step as follows

2.2.3 Extracting Naturally Arising Groups and Components

The estimated and maximised parameters can be used to guide both unsupervised and supervised modelling and in both cases, they potentially help to avoid over-fitting. We adopt PCA and data clustering for dimensionality reduction. The former seeks to transform a number of correlated variables into a smaller number of uncorrelated variables, called principal components, i.e., we can explain the variance-covariance structure of a high dimensional random vector through a few linear combinations of the original component variables. The variables in Table 2 can be formulated in a generic form as in Equation 16, where each extracted component is estimated as a weighted sum of the variables

Here, $k=1,2,3,4,5$ denotes the number of components and $i=1,2,3,4,5$ denotes the number of variables. The vectors $w_{\mathit{\text{ik}}}$ are chosen such that the following conditions are met.

1. $\Vert w_k\Vert =1$ for $k=1,2,3,4,5$

2. Each of the $\mathcal{P}{\mathcal{C}}_k$, maximises the variance $V\left\{w_k^{{}^{’}}\mathcal{P}{\mathcal{C}}_k\right\}$ and

3. The covariance $COV\left\{w_k^{{}^{’}}\mathcal{P}{\mathcal{C}}_k{w}_r{}_{}{}^{’}\mathcal{P}{\mathcal{C}}_r\right\}=0,\forall k<r$

The principal components are extracted from the linear combinations of the original variables that maximise the variance and have zero covariance with the previously extracted components. From a supervised perspective, extracted components can be viewed as new variables and from an unsupervised perspective, they may be viewed as “clusters”.

Cluster analysis is a method for grouping data according to some measures of similarity. If we assume $k$ distinct groups in our dataset, each with a specified centroid, then for each of the vectors $j=1,2,\mathrm{\dots }p,$ we can evaluate the distance from ${\mathbf{\text{v}}}_j\in \mathcal{P}{\mathcal{C}}_k$ to the nearest centroid from the set $\{{\mathbf{\text{x}}}_1,{\mathbf{\text{x}}}_2,\mathrm{\dots }{\mathbf{\text{x}}}_k\}$ as

where $d(.)$ is an adopted measure of distance, and the clustering objective will then be to minimise the sum of the distances from each of the data points in $\mathcal{P}{\mathcal{C}}_k$ to the nearest centroid. Optimal partitioning of the dataset requires identification of $k$ vectors ${\mathbf{\text{x}}}_1^{\ast },{\mathbf{\text{x}}}_2^{\ast },\mathrm{\dots },{\mathbf{\text{x}}}_k^{\ast }\in {\mathbb{R}}^n$ that solve the continuous optimisation function in Equation 18.

The solution to Equation 18 relates to the partial derivatives in Equations 7 and 8. Minimisation of the distances depends on the initial values in $\mathcal{C},$ hence if we let $z_{i=1,2,\dots ,n}$ be an indicator variable denoting group membership with unknown values, the search for the optimal solution can be through iterative smoothing of the random vector $x|(z=k),$ for which we can compute $\overline{\mu }=\mathbf{\text{E}}(x)\text{and}\delta =\left\{{\mu }_k–\overline{\mu }\right|y=k\in {\mathbf{\text{c}}}_z\}.$ In a labelled data scenario, $\{x_i,y_i\}i=1,2,\mathrm{\dots }n,$ Equation 18 transforms to the minimisation of Equation 19

where $x_i$ are described by the parameters $\left\{\overline{\mu }\text{and}\delta \right\}\in \theta$ and $g(x_i;\theta )$ are fitted values. Equations 17 to 19 relate to the K–Means clustering algorithm, which searches for clusters in numeric data based on a prespecified number of centroids. The decision on the initial number of centroids ultimately impinges on the detected clusters, and we attempt to address this issue via the two algorithms, below. Addressing spatio-temporal variations in dataset appeals naturally to dealing with randomness in data (Mwitondi and Said, 2013) and adopting interdisciplinary approaches to gain a unified understanding and interpretation of data modelling. In the next exposition, we present two algorithms: EstiMax and the Sample-Measure-Assess (SMA) algorithm (Mwitondi and Zargari, 2018), developed to address variations in data due to inherent randomness.

2.2.4 EstiMax Algorithm

Algorithm 1 below, adapted from a previous similar work (Mwitondi et al., 2018a), is designed to identify naturally arising structures involving the five pollutants in Table 2 and the associated attributes–Times, DayTimes and Period. It searches for smoothing parameters that optimise the densities in Figure 5, while minimising the effect of randomness (Mwitondi, Moustafa and Hadi, 2013).

The EstiMaxi algorithm provides general mechanics for estimating crucial parameters of the pollutants and their likelihoods. Our interest is to obtain multiple sets of parameters $\mathrm{\Theta }\{.\}$ in the sampled periods for accurate and consistent estimation. We can then associate the parameters across samples with the time-related variables in Table 2. The SMA algorithm (Mwitondi and Said, 2021) below, invokes the EstiMax, performs comparative analysis and determines the best model. However, this sequential relationship is not a requirement for implementing either of the algorithms, because the EstiMax can estimate parameters from any input dataset, while the SMA can be executed with any predetermined parameters.

The population error will always be greater than or equal to the training sample error drawn from the same population, hence the expected difference is guaranteed to be greater than or equal to 0 ($\mathbb{E}\left[\mathrm{\Delta }\right]\ge 0$). Section 3 focuses on exploratory data analyses, optimal estimation of the key parameters for the pollutants in Table 2, and their spatio-temporal associations, thus providing insights into assessing pollution patterns and relevant dynamics in southern China. The modelling techniques adopted—PCA, K-Means and the EM algorithm—are based on their underlying mechanics, particularly the influence of the starting point for the last two and their adaptability to the two algorithms.

3.2.1 Optimisation of Parameters

Figure 6 shows a 2–dimensional clustering of raw pollution datasets, which allows readers to obtain a clear visualisation of the numeric data presented in Table 2. The data points are represented by points in the ellipsoid plots, using principal components or multidimensional scaling. We use the method to illustrate the optimisation process adopted in estimating $\theta$, as defined in Equation 2. The plots represent a set of two, three and four clusters generated from average pollution levels of the five pollutants. Each of the three set of clusters is formed around a set of centroids that can be determined in a number of ways–such as random initialisation or data–dependent parameters.

Figure 6

Table 5 exhibits the contributions of each pollutant into the formation of components, as well as the final centroids for each of the three sets of clusters. In all three cases, the two components contributed 76.41% of the statistical variation.

An illustration of estimation and maximisation of key parameters of the pollutants using the EstiMax algorithm is given in Figure 7a and 7b based on raw ${\text{PM}}_{10}$ and ${\text{PM}}_{2.5}$ concentrations from the 8 stations in southern China. In each case, the dotted thin lines represent the average across samples and the thick line is the maximised estimated average across the samples. The algorithm does quite well in capturing the multi-modality of the pollutants.

Figure 7

The bimodality of the estimated pollution levels across the eight stations were the best results after comparing 2, 3, 4 and 5 modes, in relation to the raw patterns. The bimodality, effectively, implies that there are two pollution levels across the eight monitoring stations: high and low, although there are clear outliers in the upper tails in each case.

3.2.2 Extracting Components

If we denote the overall state of pollution by $\mathcal{P}$, our interest lies on the multicollinearity of the known factors, which cause inaccurate estimates of various parameters, and as a result affect future modelled outputs or forecasts. One way of addressing this issue is to apply PCA –a statistical method that reduces data dimension, and maximises variation within the same dataset. Ultimately, this results in the estimation of parameters $(\theta )$ in Algorithm 1, which improves the accuracy of $\mathcal{P}$ for being used as prior information for ongoing pollution monitoring within a prescribed spatial domain.

Figure 8a shows the correlation figures among each pair of the five pollutant variables, while Figure 8b shows the extracted components. In this case, the standard deviations in Component 1 through Component 5 are 19.26, 17.18, 13.24, 4.46 and 3.01, with proportions of variance equal to 0.426, 0.339, 0.201, 0.0228 and 0.010 respectively. Note that the first two components account for 76.5% of the variation in the pollution dataset.

Figure 8

Table 6 shows the contribution of variables in each component. It can be seen that ${\text{PM}}_{10}$ has a very strong relationship with component 4, while the absence of Nitrogen Oxide in Component 2 is most pronounced, as is that of Nitrogen Dioxide in Component 5. Similarly, the absence of ozones in Component 3 is seemingly vital in its formation.

For a clearer interpretation of the numerical and graphical results in Section 3, we need to relate them to the modelling strategy in Section 2.2–particularly to Algorithms 1 and 2. Part of the modelling strategy provides the theoretical foundations on which we can understand the overall data behaviour, and part of it provides aspects of its practical implementation. The main idea is to attain optimal interpretability of modelling results which, typically, derive from training and validating statistical models on “random samples” and ultimately applying them on new data that is also “random”, entailing variations across samples. Note that, mathematically, from an $m\times p$ dimensional dataset, a total of $p$ components can be extracted, but, usually only a few will explain the variation in the data. It is the foregoing data variation that Algorithms 1 and 2 seek to address. For instance, the centroids of the formed clusters in Table 6 are based on 1, 2, 3, 4 and 5 clusters, although for PM₁₀ and Ozone, the maximum of 75 and 114 respectively could have been extracted. Drawing new samples of PM₁₀ and Ozone from similar data sources isn’t necessarily going to generate the same number of components or the same loadings. The five clusters in Table 6 are judged optimal according to the repeated sampling via Algorithms 1 and 2. The same applies to correspondence analysis.

3.2.3 Correspondence Analysis

The average pollution levels for each of the five pollutants from each of the eight locations are given as continuous data. Their density plots exhibit multi-modality –a distributional behaviour that is supportive of discretisation of each variable. If we denote the average vector by ${\mathcal{V}}_a(x),$ we can discretise it by setting the following rule:

Equation 20 describes how an ordered continuous data vector can be visualised based on its overall behaviour. Mean estimates are in accordance with the EstiMax algorithm, and initial points can range from completely random to data-dependent, such as those in Figure 3 or other data parameters, such as using percentiles or quartiles. In each case, a quantitative vector is broken down into discrete segments that form categories of visual analysis based on correspondence analysis. The technique enables the visualisation of associations between different categories of selected data attributes in a 2–dimensional space. It seeks to establish associations between some row elements and some column elements, generating orthogonal components, with maximisation of variation in the data in mind.

Figure 9a represents associations between ${\text{PM}}_{10},$ discretised averages with the Time of Day, while Figure 9b shows the same association with annual periods. In both cases, the axes measure the levels of variation in the data, with the extreme left of the horizontal axis representing the negative measure and the extreme right representing the most positive measure. Similar definitions are applicable for the south and north directions of the vertical axis. The highest variation in both panels is accounted for by the first component, at 95.7% and 92.7% respectively.

Figure 9

Tables 7 and 8 represent the frequencies of each of the three categories during the day and the year 2023, respectively. Note how the High and Medium categories dominate in each case. It is worth mentioning that this distribution is, and should always be, conditional on the distributional behaviour of the discretised variable, since the number of classes does affect how many cases fall into each class. Our interest in Figure 9 is to find out the points that contribute to the solution provided by the method which, in this case, presents associations between the pollutants ${\text{PM}}_{2.5}$ and ${\text{PM}}_{10}$ with different time periods. CA forms these patterns based on “expected values”–associating row and column labels for the disparity between “expected” and “observed” values, in order to explain the percentage of variance in the data.

Figure 9a shows the relationship between ${\text{PM}}_{10}$ and times of day. The strength of the relationship is measured by the distance between two points and the tightness of the angle: the tighter the angle, the stronger the relationship, e.g., the angles formed between “Day” and “Evening” with “High” pollution. Orthogonal angles (${90}^{\circ }$) indicate no relationship, while a ${180}^{\circ }$ angle indicates negative association. The length of the line connecting the row label to the origin indicates the strength of the row label association, e.g., “Morning” and “Night” with respect to “Medium” and “Low” pollution levels. Being farther away from the origin means more closely associated with the factors in the proximity. In Figure 9a and 9b, most time periods are farther from the origin, without being very close to any factor.

Table 9 presents the row points, most associated with the first principal dimensions (PD) for the association between ${\text{PM}}_{10}$ and Day Times, while Table 10 presents similar information for ${\text{PM}}_{10}$ with annual periods. Evening times and January have the lowest row points association with the second PD, whereas Night Times and October have the highest association with the second dimension. These are quite interesting patterns, worth following through as they potentially indicate commonalities between those periods which may relate to specific activities or lack of them.

A further comparative analysis from eight different stations is provided in Figure 10, in which ${\text{PM}}_{10}$ and ${\text{PM}}_{2.5}$ are compared for stations 1 and 8, positioned further apart from each other. The top two panels present ${\text{PM}}_{2.5}$ concentrations, while the bottom panel presents ${\text{PM}}_{10}$ concentrations within the same stations. In both cases, variations within the four months of 2023 are more obvious than variations within different periods of a day.

Figure 10

The discretisation in Equation 20 provides the comparative basis in line 15 of the SMA algorithm (Algorithm 2), which can be visualised via plotting. On the basis of such comparison, one can update parameters within the algorithm, and conduct corresponding assessment of $\mathbb{E}\left[\mathrm{\Delta }\right],$ to determine the difference between predicted and actual parameters.