The research related to air quality is one of the important atmospheric environmental issues these days. Among them, fine dust, which is particulate pollution, is being addressed as a very significant air pollutant, as it can cause health hazards such as respiratory diseases in humans. Recently, to overcome the limitations of surface in-situ observations in the fine dust monitoring network, satellite-based measurements of aerosol optical depth (AOD) are used to estimate PM_2.5 concentrations across the world (Guo et al., 2009; Xin et al., 2014). AOD is a measure of the extinction of the solar beam by aerosols, which are tiny solid or liquid particles suspended in the air or as a gas, and measured by the Moderate Resolution Imaging Spectroradiometer (MODIS). AOD has its values typically range from 0 to 5, with higher values indicating a higher concentration of aerosol particles in the atmosphere or more effective scattering or absorption of radiation.

Various statistical models for AOD data analysis have been developed. Lee et al. (2012) developed spatial clustering models that combine AOD and ground monitoring data to predict PM_2.5 concentrations, and Karimian et al. (2023) proposed wavelet based model to forecast PM_2.5 concentrations using Himawari AOD. Various dimension reduction methods are also commonly used in analyzing AOD data. For example, Khoir et al. (2022) used rotated empirical orthogonal function (EOF) to analyze spatial and temporal AOD distribution, and Ma et al. (2021) compared EOF and nonnegative matrix factorization (NMF) in analyzing AOD over China.

In line with this latest research trend, we apply three dimension reduction methods to analyze AOD over East Asia; NMF, EOF, and ICA. Furthermore, we propose a simple, but powerful prediction method based on these dimension reduction methods.

The remainder of the paper is organized as follows: In Section 2, we describe the AOD data over East Asia. Then, we briefly review the three dimension reduction methods, NMF, EOF, and ICA in Section 3, and also propose a forecasting method. In Section 4, AOD analyzing results including prediction performances are given. Lastly, concluding remarks are presented in Section 5.

In this study, we used the MYD04 L2 dataset of Aqua MODIS (Levy and Hsu, 2015), which has a spatial resolution of 10 km and a temporal resolution of 5 minutes. We focus on the region from 110°E to 145°E in longitude, and 30°N to 45°N in latitude which covers East Asia including Korea, Japan, and eastern China (Figure 1). The number of grid points in this area is 53001. We use monthly averaged data from 2018 to 2022, and the number of time points is 60.

The left panel of Figure 1 shows the average AOD from 2018 to 2022. The highest AOD values are found in the east coastal industrialized regions of China, the western coast of Korea, and the Yellow Sea. It is well known that the winds blowing in the southwest direction have the large impacts on South Korea’s ambient air pollution levels, which is consistent with the direction of emissions from Shanghai, China (Kim, 2019). The right four panels in Figure 1 shows the average AOD in each season. We observe that the AOD value is the highest in spring (March-April-May; MAM), followed by summer (July-June-August; JJA), winter (December-January-February; DJF), and autumn (September-October-November; SON). This is consistent with the results of a study on the characteristics of AOD distribution according to the season caused by the inflow of dust from China and Mongolia and the seasonal wind (Kim et al., 2010).

3.1. Dimension reduction methods

3.2. Prediction based on dimension reduction methods

In AOD data analysis, we can obtain components related to time and spatial pattern from each dimension reduction methods. For example, in the case of NMF, W and H are time component and spatial component, respectively. That is, each row of spatial component is a latent variable made up of a linear combination for each grid of data, and each column of time component is the weight of latent variables over time. Therefore, if we can predict time component for the future time point, then AOD value for each grid point can be predicted by simply multiplying predicted time components and the corresponding spatial components. To take the case of NMF as an example, the procedure can be summurized by following algorithm.

• Given T × m dimensional data X, perform nsNMF and obtain Ŵ and Ĥ.

• For each k^th column of Ŵ, ŵ_k = (ŵ_1k, . . ., ŵ_Tk)’, predict w_T+1,k using ARIMA. Then, we have Ŵ_T+1,*: = (ŵ_T+1,1, . . ., ŵ_T+1,r) ∈ ℝ^r.

• Then, we can predict X at time T + 1 as

where j^th component of X^T+1,*NMFis computed as x^T+1,jNMF:∑k=1rw^T+1,kh^k,j for j = 1, . . ., J.

In general, h-ahead prediction also can be obtained as X^T+h,*NMF:=W^T+h,*H^, where Ŵ_T+h,* is h-ahead predicted value using ARIMA model.

Similar to NMF, we can make predictions for EOF and ICA as follows;

where V̂_{T+1, *} and Â_T+1,* are predicted coefficient vector and mixing matrix at time T+1 using ARIMA, and Ẑ and ŝ are estimated EOF matrix and source signal, respectively. Here, 1 ∈ ℝⁿ is a vector of ones and X̄ = (χ̄₁, . . ., χ̄_m)’ ∈ ℝ^m is the average value at the each grid over the training time. In EOF and ICA, we subtract X̄ from the data as a pre-processing step to make the mean of the input data equal to zero. Therefore we should add the X̄ in the prediction formula in (3.2).

Figure 2 shows the flow chart of the prediction algorithms.

4.1. Feature extraction

We apply NMF, EOF, and ICA to the AOD data, and extract features. Since there is no unified method to determine the number of components in NMF and ICA, we exam various dimension reduction level, r. As a result, r = 6 provides most interpretable results for all methods.

Figure 3 shows 6 factors of NMF and corresponding time series. We observe that the factor 2 represents the Pacific Ocean adjacent to East Asian countries and Eastern coast in China, and from ŵ₂ in Figure 3(b), the corresponding time series has seasonality peak every January. We also observe that the factor 4 represents AOD distribution of Korea, Japan and Shandong Bandao of China, and corresponding time series has seasonality peak every spring (MAM).

Figure 4 shows that 6 EOFs and corresponding PC time series, which have the largest explained variance. The 6 components explain 40.2% of the total variance. From the figure, we observe that the EOF 1 represents the East coastal industrialized regions of China and the West coastal regions of Korea, which is similar to factor 4 in NMF. From the PC 1 in Figure 4(b), the time series of factor 1 has semi-annual cycle with peaks occuring in both spring (MAM) and autumn (SON). EOF 2 represents area around Beijing and Bohai Sea, and the pattern has seasonality peak every winter (DJF).

Figure 5 shows that 6 sources and corresponding IC time series obtained from ICA. Source 1 represents the East coastal industrialized regions of China and the West coastal regions of Korea, which is similar to EOF 1 of EOF result, and factor 4 of NMF result. Source 3 represents area around Beijing and Bohai Sea, which is similar to EOF 2 in EOF result. And it has seasonality peak every winter (DJF), which is similar pattern with the PC 2 in EOF.

From the results, we observe that some components of three methods give similar spatial patterns of AOD. For example, factor 4 in NMF, EOF 1 in EOF, and source 1 all provide high AOD levels in East coastal industrialized regions of China and the West coastal regions of Korea. However, since NMF provides nonnegative coefficient matrix as shown in Figure 3(a), the results of NMF are easier to understand, as negative values hold less meaning in AOD data analysis. Also its seasonality is more clearly observed in Ŵ₁ (Figure 3(b)).

4.2. Prediction

Here, we predict monthly AOD in 2022 based on NMF, EOF, and ICA as described in Figure 2. To make prediction at time point T + 1, we use AOD data up to time point T. For example, to forecast AOD on January 2022, we use data from January 2018 to 2021 December. Then, we use data from January 2018 to 2022 January to forecast on February 2022, and so on.

The results are presented in Figure 6. We plot the average difference between the true AOD values and prediction values, ∑_t([X]_{t j} – [X̂]_{t j}), for each grid j, where [X]_{t j} and [X̂]_{t j} are the true and predicted AOD values at time t on grid j, respectively. We can observe that the difference between the true values and the all predictions are close to zero, and all three methods tend to overestimate the same region, China’s coastline. We further plot the differences for each season; spring (MAM), summer (JJA), autumn (SON), and winter (DJF) (Figure 7). For all seasons, EOF based prediction works best, but NMF based prediction works similarly except winter. On the other hand, ICA based prediction works poorly compare other two methods except winter.

For the numerical comparison, we compute root mean square error (RMSE) and mean absolute percentage error (MAPE) for each method as follows.

The results are summurized in Table 1. Here, we add simple ARIMA forecasting model for the comparison method. We apply ARIMA model for each grid, and the optimal parameters for the ARIMA model are chosen based on the Akaike information criterion (Hyndman and Khandakar, 2008). We also apply simple LSTM forecasting model as the comparison method, which is known to be good in particulate matter concentration prediction (Cho et al., 2021). The model consists of two LSTM layers in addition to one dense layer with ReLu activation function. And the model uses stochastic gradient descent (SGD) optimizer algorithm and root mean square error (RMSE) loss function.

From the results, EOF provides best performance for all seasons. ARIMA shows higher performance than NMF and ICA and lower performance than EOF for all seasons. By MAPE, EOF shows best performance in spring, summer, autumn and total period, and ICA shows best performance in winter. LSTM shows poor performance for all seasons, sine we do not have enough sample size for the learning model.

As mentioned earlier, dimension reduction level is chosen as r = 6 for the interpretation. Therefore, we further examine the result under various dimension reduction levels, r = 4, 5, . . ., 20. From the Figure 8, we observe that EOF-based method works best for all r levels. Remarkable things are that large r does not guarantee high predictive performance. For NMF case, prediction performance even gets worse as r increased, and for the case of EOF, prediction performance does not change a lot depending on r. Main reason is that first 6 components explains most significant patterns of the AOD data, and if we consider large r value, meaningless patterns are included.

We also consider 2 and 3 month ahead predictions, and Figure 9 shows the RMSE values under various r. In both 2 and 3 month ahead predictions, proposed EOF forecasting shows best performances.

In this study, we consider three dimension reduction methodologies to analyze AOD data, and compare the results in terms of feature extraction and prediction. Regional information that contributes to AOD are clearly reflected from all three methods. For example, the East coastal industrialized regions of China and the West coastal regions of Korea, which has significant contribution of AOD are well represented by all three methods. However, there is a factor in NMF that represents the Pacific Ocean adjacent to East Asian countries and Eastern coast in China, but other methods can not detect it.

Also, we predict monthly AOD in 2022 based on three methods, and compare their predictive performance. From the results, we observe that the EOF shows superior predictive performance compare to that of NMF, ICA, ARIMA and LSTM model. For various dimension reduction levels, r, EOF based method still works better than other prediction methods.

There are some limitations in this study. First, the dimension reduction level, r, has different meaning in each method. In EOF, r means the number of EOFs that explains variance of the data most, whereas r means the target number of dimension reduction level in NMF. In ICA, it means the number of statistically independent components. Various methods have been developed in determining r in each method, and we expect that we can improve the results by finding optimal number r. Second, we use simple ARIMA and LSTM model in predicting time series components in each method, but more advanced time series forecasting methods can be used in the procedure. Although these limitations, we observe that all three dimension reduction methods provide meaningful patterns of the AOD data in East Asia region, and the methods provide promising prediction performances. From the prediction performance, we conclude that the EOF has superior predictive performance over other methods.