Embedded Generative Air Pollution Model with Variational Autoencoder and Environmental Factor Effect in Ulaanbaatar City

Abstract

Air pollution is one of the most pressing modern-day issues in cities around the world. However, most cities have adopted air quality measurement devices that only measure the past pollution levels without paying attention to the influencing factors. To obtain preliminary pollution information with regard to environmental factors, we developed a variational autoencoder and feedforward neural network-based embedded generative model to examine the relationship between air quality and the effects of environmental factors. In the model, actual $S{O}_2$, $N{O}_2$, $P{M}_{2.5}$, $P{M}_{10}$, and $CO$ measurements from 2016 to 2020 were used, which were assembled from 15 differently located ground monitoring stations in Ulaanbaatar city. A wide range of weather and fuel measurements were used as the data for the influencing factors, and were collected over the same period as the air pollution data were recorded. The prediction results concerned all measurement stations, and the results were visualized as a spatial–temporal distribution of pollution and the performance of individual stations. A cross-validated $R^2$ was used to estimate the entire pollution distribution through the regions as $S{O}_2$: 0.81, $P{M}_{2.5}$: 0.76, $P{M}_{10}$: 0.89, and $CO$: 0.83. Pearson’s chi-squared tests were used for assessing each measurement station, and the contingency tables represent a high correlation between the actual and model results. The model can be applied to perform specific analysis of the interdependencies between pollution and environmental factors, and the performance of the model improves with long-range data.

1. Introduction

Serious environmental issues are common in densely populated large cities in developing countries. There are many factors contributing to this, starting with climate change; however, urbanization is a global problem [1]. In the past two decades, the number of inhabitants in Ulaanbaatar city has increased more than three times. The capital city of Mongolia sits on the Tuul River, the third-longest river in Mongolia, in a valley surrounded by large mountains, the geographical map of which is shown in Figure 1a. As the city was built in large river watershed lowland between mountains, the temperature is often inverse during the wintertime [2]. Moreover, the city’s geographical and climate features, as well as its poor heating system infrastructure, make the city one of the most air-polluted cities in the world today. The standard sewerage and pure water supply systems were initially designed for 400,000 inhabitants 46 years ago, but now around one 1.4 million people live in Ulaanbaatar, [3]. Some of these 1.4 million citizens are living in the boundary areas of the city, called ger districts in some reports, which are suburbs that are not connected to the main heating supply. In these suburbs, every family has its stoves and brick stoves as their heat source. Figure 1 shows two neighboring districts: one comprises standard heating-supplied buildings, and the other is a suburb with traditional heating stoves. Unfortunately, the suburban area increases every year [4]. Figure 1 also illustrates the suburban area and the area of standard heated buildings, with the ger district area being 6–7 times larger than the area with apartments. In the warm season, wood is the most commonly used material in traditional stoves, but in the cold season, coal is used. Furthermore, in the outer boundary areas of the ger districts, other materials are sometimes used to reach the minimum warm temperature for living, such as tires, bones, old clothes, etc.

Hence, the smoke from coal and other materials resulted in a horrifying situation in winter, as reported in some earlier works [5,6]. Before the government standardized fuel in Ulaanbaatar, in one example, the harmful dust level was 6–7 times higher than the most lenient World Health Organization standards, according to the World Bank Report in 2012 [7].

From 15 May 2018, the Government of Mongolia strictly banned the import and use of raw coal in Ulaanbaatar, and issued an order for the use of improved fuel instead. Hence, the air-quality index (AQI) has improved, and there was a big difference in the AQI between 2019 and 2020, especially in the wintertime. Consequently, the AQI is influential in addressing air-quality improvement because climatology contributed to improving fuels. However, the role of the AQI measurement network needs to be clarified. To be precise, even devices that accurately measure air quality are only informed by values of measurement from the past, with very few cases in real-time. Alternatively, the tools are simply about listening, i.e., how harmful the AQI was. In this situation, a predictive method may be the solution regarding the historical data information of the AQI, climatology, and improved fuel. Thus, we developed an embedded generative model for the AQI that depends on the environmental factor effect (EFE) influencing the AQI, which is a result of the city’s geographical features, such as temperature, wind, atmospheric pressure, etc., including improved and raw fuel measurement quantities.

Therefore, studying the interdependence between AQI and EFE with fuel measurement quantities is the purpose of this work. The variational autoencoder (VAE) and multilayer perceptron methods were used to develop an embedded generative approach. From the VAE learning, we produced new data based on its latent variables. Latent variables play a crucial role in VAE generation, highlighting data and reducing its dimensionality in a shallow space. The EFE with fuel measurement quantities was measured simultaneously with the AQI. The dimensions of the latent variables and the EFE with fuel measurement quantities were the same size. Hence, to discover the relationship between the AQI and the EFE data’s latent variable, a multilayer perceptron was used. Thirty-two variables were gathered in the EFE data, such as ambient temperature, pressure, and wind speed, as well as the impact of raw coal and improved fuel information. For the VAE learning process, a four-year period from 11 of Ulaanbaatar’s (UB’s) air-quality measurement networks were used. The proposed method’s results were evaluated by a Pearson chi-squared test and compared with a non-linear regression prediction in the final step.

The measurement quantities of raw coal and improved fuels in three different laboratories are shown in

Table 1. Fuel quality comparison.

Fuel Type	Measured Laboratory	Q1	Q2	Q3	Q4	Q5	Q6	Q7	Q8	Q9
Improved Fuel	Southwest region, Russia	2.4	22.2	19.6	0.89	6374	67.63	3.81	1.83	3.68
	Central Geological Laboratory	1.9	22.7	19.8	0.87	6984	-	-	-	-
	Laboratory of Mineral Resources and Petroleum Authority	2.7	23.3	19.3	0.86	6334	65.1	5.4	1.32	7.57
	Institute of Chemistry and Chemical Technology	0.8	22.9	18.68	0.89	5918	63.86	3.76	1.69	-
	National standard MNS 5679:2019	≤10	≤29	≤22	≤1.0	≤4200	-	-	-	-
Raw Coal	National standard MNS 3818:2011	37.5	17.5	44.8	0.38	3360	-	-	-	-
	National standard MNS 6226:2011	6	28	26	1.3	5500	-	-	-	-
	CRRI China	-	-	-	-	-	67.3	4.2	0.95	17.2
	CRRI China	-	-	-	-	-	69.3	3.8	3.8	3.1

. The items in Table 1 are as follows: Q1—received basis refers to the moisture content of the coal (maximum percentage of water); Q2—dry basis ash content (maximum percentage of dry); Q3—dry basis volatile matter (maximum percentage of dry); Q4—total dry basis sulfur content (maximum percentage of dry); Q5—lower caloric value (kcal/kg min); Q6—percentage of carbon; Q7—percentage of hydrogen; Q8—percentage of nitrogen; and Q9—percentage of oxygen.

Related Work on ANN and Deep Learning for Ambient Air Pollution Estimation

The application of a neural network (NN) for air pollution is associated with the development of an artificial neural network (ANN), with mainly machine learning and VAE being applied in prediction works. For instance, in early works, the ANN technique was applied for pollution prediction: Reich et al. used 24 h $S{O}_2$ data for simple ANN architecture and discussed some of the limitations of the approach [8]. Moustris et al. used much larger meteorological and air pollution data to predict the pollution index of $N{O}_2$, $CO$, $S{O}_2$, and $O_3$ three days ahead in seven different places in Greece [9]. For urban air pollution prediction, an ANN was applied and compared with multiple time series methods, such as ARIMA, fuzzy time series, and principal component analysis (PCA), using 7–10 years of Malaysian air pollution measurement data [10,11,12]. A more expanded description of the ANN-based approach can be found in the review work in [13,14,15]. Furthermore, an ANN was used for the integrated assessment modeling of global climate change and ecological research [16,17].

Machine learning (ML) approaches have been applied for more automatic and accurate air pollution predictions from a large amount of input data with numerous outputs. Zhang et al. studied the superior predictive ability of ML methods based on six years of Hong Kong air-quality index measurements [18]. Moreover, the special gradient-boosting ML approach was applied for predicting the $P{M}_{2.5}$ concentration in China [19]. Delavar et al. used an autoregressive non-linear neural network with an external input for improving air pollution prediction in the case of Tehran city [20]. The research of Ly et al. considered the dependence between the pollution index and the sources obtained from multisensory measurements and weather data in air pollution predictive analysis-based ML approaches [21]. Similar work was reported by Arnaudo et al. in the urban area of Milan city using different configurations employing machine and deep learning models, such as a linear regressor, an artificial neural network using Bayesian regularization, a random forest regressor, and a long short-term memory network [22].

Lee et al. presented a gradient boosting–based machine learning approach for predicting $P{M}_{2.5}$ concentration based on a large-scale data set in Taiwan city [23]. For UB city, $P{M}_{2.5}$ concentration prediction has been studied using two different ML approaches, and high-performance results have shown the spatiotemporal variations in the high-emission areas [24]. Moreover, Meredith et al. compared several machine learning methods for ground-level pollutant mixtures of particulate matter with a diameter smaller than $P{M}_{2.5}$ for air-quality assessment in Mongolia [25]. The survey works of ML methods for ambient air pollution assessment are summarized in [26,27,28,29,30].

Among the generative models for air contamination, deep learning methods are a more recently developed approach. Tien et al. worked on air-quality inference based on a graph-based matrix completion problem, and applied a variational model on graph convolutional autoencoders [31]. Additionally, they used a mobile Internet of Things sensor for fine-grained air quality in deep learning inference [32]. The spatiotemporal prediction of air pollution based on deep learning approaches has been studied widely for many different purposes. Li et al. studied the spatiotemporal deep learning–based air-quality prediction method that inherently considers the spatial and temporal correlations and prediction results of all stations compared with traditional time series prediction models [33].

Li et al. introduced the advantage of using two-stage models in geographically weighted machine learning [34], and deep learning to robustly impute a long time series of multi-angle implementation of atmospheric correction aerosol optical depth [35]. The achievements of embedded methods with deep learning approaches were considered in [36,37], and the high performance of autoencoder-based residual networks consisting of learning models is analyzed in [38].

Fan et al. studied the framework and feasibility of an idea based on deep RNN [39], whereas improved feature analysis was considered in [40] based on fine-grained air quality.

2. Methodology

2.1. Area of Study and Data Collection

The capital city UB has a history spanning more than 360 years and is currently home to 50% of the country’s population. This population is continuously growing; therefore, the most significant problem is in an urban ecological context, with environmental pollution also dramatically increasing. In addition, the city has a harsh and unique climate, and its geographical location behavior is the main cause of this air pollution. The city is located in a large valley surrounded by mountains, at an elevation of approximately 1580 m above sea level. Due to the elevation and being located far from the sea, the city has a very long-lasting winter, making UB one of the coldest capital cities in the world. Ulaanbaatar receives an average of 239 mm of precipitation per year, of which 1 mm or more falls in 39 days, and 5 mm or more of rain falls in 15 days. Generally, this area is arid, with the humidity sometimes decreasing to 4% per year. The wind comes mainly along the valley, but the strength of the wind is low throughout the year, with 10% of the wind observed in the spring months.

The average annual surface air temperature is −0.7 °C, dropping to −45 °C in the winter and warming up to +59 °C in the summer. Sudden frosts occur on the surface of the soil in all months except July. Therefore, the ground freezes to a considerable depth of 4–8 m. The average wind speed for the country is less than 2 m/s, which affects the air pollution levels, making them more constant for an extended time. Additionally, the average wind speed of the capital city is primarily 0.5 m/s from the northern and northwestern areas, which means that there is not much air circulation. Thus, air pollution remains stable for long periods around the center of the city. Figure 2 shows the wind information maps across four years.

The National Agency for Meteorology and Environmental Monitoring (NAMEM) established the first air-quality monitoring station in Ulaanbaatar in 1977. This was the foundation for the establishment of an air-quality monitoring network in Mongolia. Until 2008, there were four permanent air-quality monitoring stations in Ulaanbaatar, which measured only two indicators, sulfur dioxide and nitrogen dioxide. There was no continuous air-monitoring system built on the outskirts of Ulaanbaatar city.

Due to Ulaanbaatar’s air pollution issues, there is a need to expand the network by equipping it with automatic, continuous measurement capabilities and improving its capacity. The air-quality monitoring network consists of the Air Pollution Reducing Department of Capital City (APRD) Network and the NAMEM Network. A German grant established the APRD Network with five stations in 2008, and the NAMEM Network, with 10 stations, was established by means of a French loan in 2010. Currently, 12 automatic stations are measuring sulfur dioxide, nitrogen oxides, carbon monoxide, ozone, $P{M}_{10}$, and $P{M}_{2.5}$ every 15–30 min, and 3 non-automatic controls of sulfur dioxide and nitrogen dioxide (Figure 3). The concentration is determined using a chemical solution.

2.2. VAE for Embedded Generative Method

Let $\mathit{y}=\left\{y_1,y_2\dots y_n\right\}$ be the AQI measurements collected from measurement networks at different locations in the city. Furthermore, $\mathit{x}=\left\{x_1,x_2\dots x_n\right\}$ represents the EFE of the air quality, which is measured at the same time as $\mathit{y}$. Next, the purpose is to understand the relationship between those two quantities; $F:\mathit{x}\to \mathit{y}$. In the case of air quality, the relationship is completely unknown; generally, it is assumed to be a system of highly non-linear maps. Moreover, data of $\mathit{x}$ and $\mathit{y}$ raw and sampling methods are very different in each measurement. In this situation, a variational autoencoder (VAE) has the benefit of representing the data in low dimensional continuous latent variables with intractable posterior distributions [41].

The VAE method optimizes the weight parameters $\mathbb{W}$ and bias $\mathrm{𝕓}$ of the function: $\overline{F}\left(\mathbb{W},\mathrm{𝕓},x\right)$, which is a non-linear approximation of the relation map, as follows:

$$F=\underset{\mathbb{W},\mathrm{𝕓}}{\mathrm{argmin}}\frac{1}{N}{\displaystyle \sum }_{i=1}^N\parallel {\mathit{y}}_{\mathit{i}}-\overline{F}\left(\mathbb{W},\mathrm{𝕓},{\mathit{x}}_{\mathit{i}}\right){)\parallel }^2$$

(1)

where ${\left\{{\mathit{y}}_{\mathit{i}},{\mathit{x}}_{\mathit{i}}\right\}}_{i=1}^N$ is the data set provided and stored from the relation map $F$. The optimization process of $\overline{F}\left(\mathbb{W},\mathrm{𝕓},{\mathit{x}}_{\mathit{i}}\right)$ consists of two maps, encoder Φ and decoder Ψ. The encoder reduces the dimension of y into dimensionally reduced space, Φ: $\mathit{y}$ → $\mathit{h}$, and then the decoder map generates new $\widehat{\mathit{y}}$ values from the reduced space, Ψ: $\mathit{h}$ → $\widehat{\mathit{y}}$. Thus, when an optimal parameter is reached, we obtain the pair of $\left(\Phi , \Psi \right),$ which satisfies the following minimization:

$$\left(\Phi , \Psi \right)=\underset{\mathbb{W},\mathrm{𝕓}}{\mathrm{argmin}}\frac{1}{N}\left[ {\displaystyle \sum }_{i=1}^N\parallel {\mathit{y}}_{\mathit{i}}-\Psi \left(\Phi \left({\mathit{y}}_{\mathit{i}}\right)\right){)\parallel }^2+KL\left(N\left(\mu ,\sigma \right),N\left(0,I\right)\right)\right]$$

(2)

where $KL$ is the Kullback–Leibler (KL) divergence loss between $N\left(\mu ,\sigma \right)$ and $N\left(0,I\right)$(along with this formulation, we refer [42,43]). KL measures how the distribution of $N\left(\mu ,\sigma \right)$ associated with data $\mathit{y}$ is different to the normal distribution, defined as:

$$KL\left(N\left(\mu ,\sigma \right),N\left(0,I\right)\right)=\frac{1}{2}{\displaystyle \sum }_{s=1}^l\left({\mu }_s^2+{\sigma }_s^2-1-\ln ({\sigma }_s^2\right))$$

(3)

where l is the dimension of dimensionally reduced space $\mathit{h}$, called latent space. To obtain the latent space from the data, the multilayer perceptron (MLP) is the advised method. The map of reducing dimension consists of several hidden layers, and every layer has finite nodes. In the encoder process, the number of nodes in the steps of layer is conventional in reducing the order. $w_{ij}^k$ is denoted as the weight parameters of the layers between kth and (k + 1)th, whereas i and j are node positions, respectively. Then, the kth layer weight parameters can be represented in the following matrix form:

$$W^{\left(k\right)}=\left[\begin{array}{ccc}{w}_{1,1}^k& \cdots & w_{m,1}^k\\ ⋮& \ddots & ⋮\\ w_{1,n}^k& \cdots & w_{m,n}^k\end{array}\right]$$

(4)

where m is the number of notes in the kth layer, and n is the number of nodes in the (k + 1)th layer. To curtail the dimension of input that the number of nodes in the layer satisfies:

$$N>\cdots \underset{k-1}{\underbrace{\cdots >}}\stackrel{k}{\overbrace{m>n}}\underset{k+1}{\underbrace{>\cdots }}\cdots >l$$

In recurrent form, this can be described as:

$${\mathit{h}}^{\left(1\right)}=W^{\left(1\right)}\mathit{y}+b^{\left(1\right)}$$

(5)

$${\mathit{h}}^{\left(k\right)}=W^{\left(k\right)}{\mathit{h}}^{\left(k-1\right)}+b^{\left(k\right)}$$

(6)

where $b^{\left(k\right)}={\left[b_1^{\left(k\right)}, b_2^{\left(k\right)},\dots ,b_n^{\left(k\right)} \right]}^{⊺}$ refers to the bias parameters and $⊺$ is the matrix transpose operation. Every $h^{\left(k\right)}=\left\{h_1^{\left(k\right)}, h_2^{\left(k\right)},\dots ,h_n^{\left(k\right)}\right\}$ retains the key feature of the previous steps, and the final space $\mathit{h}$ is considered to be the feature space. Then, the encoder function is represented as:

$$\Phi \left(y\right)=\varsigma \left(W^{\left(k\right)}\dots \varsigma \left(W^{\left(2\right)}\varsigma \left(W^{\left(1\right)}y+b^{\left(1\right)}\right)+b^{\left(2\right)}\right)+\dots +b^{\left(k\right)}\right)$$

(7)

where $\varsigma$ is the sigmoid function. The main feature that distinguishes a VAE from other generating models is its latent space. If we look at latent space alone, its actual value is a quantity that cannot be measured directly. However, it contains the core of the essential information in the data, determining in the latent space the KL-divergence loss and reparameterization trick is vital. In the process of optimization, Equation (3) provides a close distribution to normal with $\mu$ and $\sigma$ expressed as:

$$\mu =W^{\mu }\Phi \left(\mathit{y}\right)+b^{\mu }$$

(8)

$$\sigma =W^{\sigma }\Phi \left(\mathit{y}\right)+b^{\sigma }$$

(9)

where $W^{\mu }$, $W^{\sigma }$, $b^{\mu }$, and $b^{\sigma }$ are the weight and bias for $\mu$ and $\sigma$. Then, the latent variable is described as:

$$\mathit{h}=\mu +\sigma ⨀ϵ$$

(10)

where $ϵ$ is the normally distributed auxiliary noise variable and $⨀$ is the element-wise product. Figure 4 shows the reparameterization trick on $\mathit{h}$, the difference with and without $ϵ$ is intractable posterior distributions.

More precisely, the latent space is centered at the coordinate and distributed in a bounded interval in any direction; thus, any linear combination between two points of the space belongs to its space (for more detailed properties of this phenomena, refer to Section 3.4 of [42]).

When $\mathit{h}$ variables are found, the decoder function is defined as:

$$\Psi \left(h\right)=\varsigma \left(W^{2k}\dots \varsigma \left(W^{k+2}\varsigma \left(W^{\left(k+1\right)}h+b^{\left(k+1\right)}\right)+b^{\left(k+2\right)}\right)+\dots +b^{\left(2k\right)}\right)$$

(11)

Then, the entire process of VAE is to determine the parameters of $W$ and $b$, where:

$$\mathbb{W}=\left\{W^{\left(1\right)},\dots ,W^{\left(k\right)},W^{\left(\mu \right)},W^{\left(\sigma \right)},W^{\left(k+1\right)},\dots ,W^{\left(2k\right)}\right\}$$

(12)

$$\mathrm{𝕓}=\left\{b^{\left(1\right)},\dots ,b^{\left(k\right)},b^{\left(\mu \right)},b^{\left(\sigma \right)},b^{\left(k+1\right)},..,b^{2k}\right\}$$

(13)

The process of optimizing those parameters is based on backpropagation and a gradient descent algorithm. Initially, it starts with random choice in $\mathbb{W}$ and $\mathrm{𝕓}$, then the sensitiveness of training errors in little changes in those parameters, and the derivative of the training error leads to the chain rule and backpropagation. The network architecture is represented in Figure 5, and the architecture of the VAE, and the algorithm of optimizing the parameters, are described in the next section, along with the embedded network.

2.3. Embedding Network for the Relationship between EFE and AQI

The formulas described in Equations (2)–(13) are basic VAE methods for generating new AQI values from the latent variables obtained in the learning process. In order to approximate the function F, whereby $\mathbb{W},\mathrm{𝕓}$, and $\mathit{x}$ are as described in Equation (1), we need different approaches to the typical VAE approach. However, the method defined in Equations (2)–(13) determines the vital space of the AQI, which is the latent variable $\mathit{h}$. To find the relationship between latent space data and the EFE, we use the feedforward neural network approach, as the period of sampling time is equal.

Expressly, the first column of Equation (14) represents n daily measurements of AQI, and the next column is the latent variables obtained through the algorithm applied to daily measurements of the AQI. Additionally, EFEs were also measured on the same day with $\mathit{y}$:

$$\begin{array}{ccc} {\mathit{y}}_1=\left\{y_1^{\left(1\right)},\dots ,y_m^{\left(1\right)}\right\}& {\mathit{h}}_1=\left\{h_1^{\left(1\right)},\dots ,h_l^{\left(1\right)}\right\}& {\mathit{x}}_1=\left\{x_1^{\left(1\right)},\dots ,x_l^{\left(1\right)}\right\}\\ {\mathit{y}}_2=\left\{y_1^{\left(2\right)},\dots ,y_m^{\left(2\right)}\right\}& {\mathit{h}}_2=\left\{h_1^{\left(2\right)},\dots ,h_l^{\left(2\right)}\right\}& {\mathit{x}}_2=\left\{x_1^{\left(2\right)},\dots ,x_l^{\left(2\right)}\right\}\\ \dots \dots \dots & \dots \dots \dots & \dots \dots \dots \\ \underset{\mathit{y}}{\underbrace{{\mathit{y}}_n=\left\{y_1^{\left(n\right)},\dots ,y_m^{\left(n\right)}\right\}}}& \underset{\mathit{h}}{\underbrace{{\mathit{h}}_n=\left\{h_1^{\left(n\right)},\dots ,h_l^{\left(n\right)}\right\}}}& \underset{\mathit{x}}{\underbrace{{\mathit{x}}_n=\left\{x_1^{\left(n\right)},\dots ,x_l^{\left(n\right)}\right\}}}\end{array}$$

(14)

Hence, the data are available to learn the relationship between $\mathit{h}$ and $\mathit{x}$. Apply the feedforward neural network to approximate the function $f:\mathit{h}\to \mathit{x}$; indeed, it is the universal approximation theorem for $f$. Thus, it gives the $f\left(\mathit{𝔀},\mathit{𝓫}; \mathit{x}\right)$ function, which is an approximation function as $f\sim f\left(\mathit{𝔀},\mathit{𝓫}; \mathit{x}\right)$ with weight and bias parameters corresponding to network nodes. Then, the new generating AQI values can be described through the embedding as:

$$\{\begin{array}{c}f:\mathit{x}\to \mathit{h}\\ \Psi :\mathit{h}\to \widehat{\mathit{y}}\end{array}$$

(15)

The architecture of the two combined networks is represented between $\mathit{y}$, $\mathit{h}$, and $\mathit{x}$ in Figure 6.

The Algorithm 1 for the embedded generative model:

Algorithm 1: embedded modeling algorithm

Result: find the relation map ƒ; ($\mathit{𝔀}$,$\mathit{𝓫}$)-initialization; ζ-step size hyperparameter; set error estimation; set algorithm for ($\mathit{𝔀}$,$\mathit{𝓫}$) same with ($\mathbb{W}$,$\mathrm{𝕓}$);

Result: embedding of Ψ, ƒ; Generate new AQI values based on the embedding; Ψ: (ƒ: x $\to$ h) $\to$ ŷ

3. Experiment

3.1. Data Expression

To validate the proposed method of illustrating the interdependency between AQI and EFE, we used four years of air pollution concentration measurements and four years of weather data. The correlation between AQI and pollution concentration is defined as:

$$AQI=\frac{(I_{high}-I_{low})}{(C_{high}-C_{low})}\left(C-I_{low}\right)+I_{low}$$

(16)

where $C$ is the amount of pollutant concentration per cubic meter of air, $C_{low}$ and $C_{high}$ are the lowest and highest breakpoints of $C$, and $I_{low}$ and $I_{high}$ are the index breakpoints corresponding to $C_{low}$ and $C_{high}$. However, the AQI is designed to disseminate air-quality information using color coding as it is accessible to the general public, normalizing the measurement into world standards, with the relationship between AQI and $C$ being direct one-to-one. Furthermore, the pollutant concentration is available from measurement stations, with most of them even working automatically, and we used $C$ instead of AQI. In addition, the pollution concentration data that were assembled between 1 January 2016 and 1 June 2020 from 15 pollution measurement stations and positions are shown in Figure 3. Three of them determined the pollution concentration based on active sampling tests, and the remaining stations automatically measured the concentration. The automatic stations used the light source to obtain the signal difference between the reference and the incoming air in sampling manifolds. This working principle analyzed the amount of different pollution concentrations by micrograms per cubic meter in a minute and eliminated rotation for optic adjustment. Related to the working principle and purpose, each station had different measuring capabilities, namely $S{O}_2$ measured from 15 stations, $N{O}_2$ measured from 13 stations, $P{M}_{10}$ measured from 14 stations, $P{M}_{2.5}$ measured from 8 stations, and $CO$ measured from 10 stations, with a total of 60 pollution concentrations gathered for the experiment. As in the example, the daily concentration of $S{O}_2$ in 2016 is expressed in Figure 7, which was measured by UB-2, with the tolerable content level being $50 \mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$.

The EFE data were composed of measurements of a comprehensive range of weather data sampled from five automatic weather measurement stations in UB. Weather data involve various levels of wind speed, wind direction, temperature, relative humidity, and the most crucial information for the fuel period. Furthermore, weather information was collected within the same period as the air pollution concentration measurement.

Surface wind regimes depend on local climatic conditions and atmospheric circulation. Each season, the average wind speeds are very different in UB, $1.8–6.3 \mathrm{m}/\mathrm{s}$ in spring, $-2.5 \mathrm{m}/\mathrm{s}$ in winter, and $1.2–5.4 \mathrm{m}/\mathrm{s}$ in other seasons. In winter, due to the strong anti-cyclones, the prevailing wind speed slows down, whereas in other seasons, the west, northwest, and north winds prevail along the prevailing air flow. Depending on the location of the weather stations in UB, the wind direction, speed, and frequency are quite varied. For example, in the station at the airport, the east wind is predominant in all seasons except spring, different to the station at the university. The four years of the probability distribution of the direction and wind speed are visualized in Figure 8. The degree of wind direction is defined in cardinal degrees as compass points.

Furthermore, four years of temperature and relative humidity information are expressed in Figure 9. For both quantities, the frequency is counted with respect to amounts of years, and the absolute measurements are represented in a bivariate histogram. In the illustration of temperature, the higher frequencies appear more often in the cold season over the four years of data than in the warm season. The amount of relative humidity is inversely related to air temperature and depends on cloud precipitation. The average annual relative humidity around Ulaanbaatar is 55–76%.

3.2. Visualization of the Results

In the VAE-based learning process, to obtain the feature space $\mathit{h}$ of the AQI data, we need to separate data into two parts, i.e., training and testing. The data size was 365 $\times$ 60 = 87,600 from 8–15 air pollution concentration measuring stations, and the data were separated into 60,000 and 27,600 for training and testing, respectively. The computation architecture for the VAE consisted of three hidden layers with nodes 40, 30, and 20, and one input channel had a length of 60. The training epoch was 200,000, the batch size was 1000, and the learning rate was 0.001. We used an RmsProp optimizer with the ReLU activation function in the hidden layer and the sigmoid activation function in the output layer.

Python and TensorFlow were used for the VAE implementation, and a PC equipped with an Intel(R) Core (TM) [email protected], with an NVIDIA GeForce GTX 1060ti 6GB GPU was used for the computation. The Python plotting tool, MATLAB, and LaTeX pgfplot packages were used in the data and result visualization. In Figure A1, Figure A2 and Figure A3 in the Appendix A, the typical air pollution measurement of $S{O}_2$ and the convergence rate of the VAE computation at a particular epoch are shown.

With this latent space computation, we completed the VAE using the algorithm illustrated in pseudocode 1. Next, for the relationship between the EFE data and the latent variables of the VAE, the development of a feedforward NN for learning was required, which is observed in the second step of pseudocode 1. The daily measurement of the EFEs consisted of 32 quantities. Hence, the feedforward NN’s input channel was 32, and it had one hidden layer with 25 nodes. Along with a training epoch of 10,000, a batch size of 20, and a learning rate of 0.001, we used the Adam optimizer with the ReLU activation function in the hidden layer and the sigmoid activation function in the output layer, whereas 70% and 30% of the data were used for the training and testing sets, respectively.

Finally, to validate our suggested method, the EFE data of the first three months of 2020 were selected for generating the predicted air pollution concentration. With respect to the first 90 days and the pollution measurement stations, the model provides 90 $\times$ 60 = 5400 as the number of results. This large number of results requires an expression of the outcome in an optimal way, and we elected a few specific depictions. For this reason, the most frequently polluted days of each pollution index were chosen to illustrate the results.

Table 2. The most frequently polluted air pollution dates in 2016–2019.

	$\mathit{S}{\mathit{O}}_2$	$\mathit{N}{\mathit{O}}_2$	$\mathit{P}{\mathit{M}}_{10}$	$\mathit{P}{\mathit{M}}_{2.5}$	$\mathit{C}\mathit{O}$
Jan	12	11	11	11	12
Feb	1	1	3	4	1
Mar	1	1	31	2	1

represents the most frequently polluted dates in the first three months of our four years of data.

The coldest days of the year directly affected the reason for representing the same behavior for gaseous pollution in the first three months of integration in Table 2. By long-standing average, the coldest days of the year are between 8 and 15 January. Therefore, in January, $S{O}_2$, $N{O}_2$, $P{M}_{2.5}$ and $CO$ are formed from fuels burned in home stoves for heating in the ger districts with high atmospheric pressure and temperature inversion. However, particular pollution is related to poor infrastructure. When the warm season starts, particular pollution from the ger districts increases due to the mobility of a large number of vehicles; consequently, at the end of March, the measurement of $P{M}_{10}$ greatly increases.

Furthermore, the inverse distance weight (IDW) interpolation method was used for the pollution concentration map. The inverse weight is defined as:

$${\omega }_{ij}=\frac{1/d_{ij}^2}{{{\displaystyle \sum }}_{i=1}^{n}1/d_{ij}^2}$$

(17)

where $n$ is the number of points for interpolation, and $d_{ij}$ is distance between the $i$th interpolation point and the jth air pollution measurement station. Then, based on Equations (14) and (15), the concentration values at the $i$th interpolation point is estimated as:

$${\mathit{y}}_i={\omega }_{ij}\Psi \left(f\left({\mathit{x}}_j\right)\right)$$

(18)

The results of Equation (18) are expressed in Figure 10, Figure 11 and Figure 12, corresponding to the integrated values in Table 2; 1 February and March for $S{O}_2$ and $CO$, and 11 January for $P{M}_{2.5}$ were selected, respectively.

For more results, we randomly selected one of the pollution measuring stations for every five air pollutions, and 90 days of both actual and model prediction results are compared in Appendix A. To evaluate the performance of the proposed model, the most commonly used correlation coefficient and Pearson’s chi-squared methods were used. The positive correlation coefficient indicates higher correlation between testing values. As in the example, Figure 13 and Figure 14 represent the correlation coefficients between the actual and predicted values of $S{O}_2$, $P{M}_{2.5}$, $P{M}_{10}$, and $CO$ obtained from all measurement stations, respectively. In the figure, the corresponding correlation coefficients are highlighted in red in the correlation matrix plot. For more validation purpose, we randomly selected one measurement station and computed the contingency tables for 90 days of actual and predicted values, which expressed in Figure A4, Figure A5, Figure A6, Figure A7 and Figure A8 in Appendix A. Two of the corresponding contingency tables for $N{O}_2$ and $P{M}_{2.5}$ are shown in Figure A9.

For the remaining three air pollution concentrations, the contingency table was similar, whereas the average chi-squared result was ${\chi }_{ave.}^2=115.14$. Thus, the statistics for the Pearson’s chi-squared test were highly dependent on each other.

4. Conclusions

In this study, a variational autoencoder-based embedded generative air pollution model was developed. The proposed approach represented substantial progress in the reliable prediction of five air pollution concentrations. The embedded generative model consists of a VAE and a feedforward NN to discover the relationship between the air pollution concentration and environmental factor effects. The convergence of learning networks based on the gradient descent method with RmsProp and Adam optimizer techniques were applied, and the convergence rate was $1\times {10}^{-3}$, consequentially. To construct the proposed model, we simultaneously assembled four years (2016–2020) of air pollution concentration measurements and weather data. In addition, information on improved and raw fuel was involved in the environment factor data. The prediction values of the model corresponded to all the pollution measurement stations, with a total of 15 $\times$ 90 = 1350 prediction values obtained for the test with respect to 15 stations in the first 90 days of 2020. Thus, the results can be visualized in two ways. First, a pollution map can be created by using the inverse distance weight interpolation throughout the regions of UB, which expresses a daily pollution map for estimating the pollution concentration for local coordinates of the city. For instance, the difference between the local measurement and the prediction results for $P{M}_{10}$ in ger district areas located in UB suburbs in four directions was 21.24 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$, −9.32 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$, 11.05 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$, and −11.87 $\mathsf{\mu }\mathrm{g}/{\mathrm{m}}^3$, respectively, in the Airport, Zuragt, Tolgoit, and Amgalan districts on 16 January 2020. The correlation coefficient method was applied while constructing the pollution concentration map of the city for each pollution index, and the average of these indicates the entire performance of the model as 0.8280. The Pearson’s chi-squared test method was used for testing the second visualization performance of the individual pollution measuring stations and expressed a high correlation between actual and predicted. A further aim of considering the VAE and NN embedding to investigate inversely the interdependencies between the weather and the fuel conditions, by sampling normal to good days of air pollution, could be considered for future study.