Spatial Gap-Filling of Hourly AOD Data from Himawari-8 Satellite Using DCT (Discrete Cosine Transform) and FMM (Fast Marching Method)

Abstract

Since aerosol has a relatively short duration and significant spatial variation, satellite observations become more important for the spatially and temporally continuous quantification of aerosol. However, optical remote sensing has the disadvantage that it cannot detect AOD (Aerosol Optical Depth) for the regions covered by clouds or the regions with extremely high concentrations. Such missing values can increase the data uncertainty in the analyses of the Earth's environment. This paper presents a spatial gap-filling framework using a univariate statistical method such as DCT-PLS (Discrete Cosine Transform-based Penalized Least Square Regression) and FMM (Fast Matching Method) inpainting. We conducted a feasibility test for the hourly AOD product from AHI (Advanced Himawari Imager) between January 1 and December 31, 2019, and compared the accuracy statistics of the two spatial gap-filling methods. When the null-pixel area is not very large (null-pixel ratio < 0.6), the validation statistics of DCT-PLS and FMM techniques showed high accuracy of CC=0.988 (MAE=0.020) and CC=0.980 (MAE=0.028), respectively. Together with the AI-based gap-filling method using extra explanatory variables, the DCT-PLS and FMM techniques can be tested for the low-resolution images from the AMI (Advanced Meteorological Imager) of GK2A (Geostationary Korea Multi-purpose Satellite 2A), GEMS (Geostationary Environment Monitoring Spectrometer) and GOCI2 (Geostationary Ocean Color Imager) of GK2B (Geostationary Korea Multi-purpose Satellite 2B) and the high-resolution images from the CAS500 (Compact Advanced Satellite) series soon.

1. Introduction

Aerosol is defined as the particles in solid and liquid floating in the atmosphere and affects the mechanism of radiation transfer and cloud microphysics (Twomey, 1974; Albrecht, 1989). Because it can deteriorate atmospheric visibility and human health (Wang et al., 2009; Zanobetti and Schwartz, 2009), a quantitative understanding of aerosol distribution is essential; South Korea also makes many efforts for ground observation and remote sensing of aerosol (Lee et al., 2020). Aerosol has a relatively short duration and significant spatial variation (Jin et al., 2018), so the point-based ground observations are not sufficient to understand the aerosol characteristics in a spatially continuous scale (Higurashi et al., 1999). Hence, satellite observations become more important for the spatially and temporally continuous quantification of aerosol. Long-term satellite data can be used for the analyses of climate change and air quality trends (Kafuman et al., 2005; Hoff et al., 2009; Bae et al., 2017).

AOD (Aerosol Optical Depth) or AOT (Aerosol Optical Thickness) is the quantity to represent aerosol using the integral extinction coefficient (NMSC, 2012), the degree that aerosol disperses the radiation along the path between a satellite sensor and the land surface, by absorption and scattering of the radiation (Kinne et al., 2006). Polar-orbiting satellite sensors such as MODIS (Moderate Resolution Imaging Spectrometer), VIIRS (Visible Infrared Imaging Radiometer Suite), MISR (Multi-angle Imaging Spectroradiometer) conduct AOD observations one or twice a day (Hsu et al., 2006; Levy et al., 2007; Hsu et al., 2013; Kalashnikova et al., 2013; Levy et al., 2013; Zhang et al., 2016; Garay et al., 2017), but tracking of diurnal variation of AOD is difficult (Gao et al., 2021). Meanwhile, recent geostationary satellites such as the American GOES-17 (Geostationary Operational Environmental Satellite) (Schmit et al., 2005), the Japanese Himawari-8 (Bessho et al., 2016), and the Korean GeoKompsat-2 (Jee et al., 2020) provide hourly AOD products during daytime. They have more advanced spatial and temporal resolutions, spectral bands, and SNR (Signal to Noise Ratio) (Schmit et al., 2017). The performance of the retrieval algorithm was improved for a spatially and temporally continuous monitoring of AOD at an interval of 10 minutes on a 2 km grid (Lee et al., 2020).

South Korea, which is affected by the westerlies, may have a significant aerosol inflow from the continent. Citizen’s worries are increasing because of the recent serious concentration of particulate matter (Kim et al., 2021). More than 300 Air Korea stations have been set up to provide point-based real-time information on air quality nationwide. However, spatially continuous understanding of the transport of particulate matter from the continent to the Korean peninsula is still difficult, so the studies of monitoring with CTM (Chemical Transport Model) and remote sensing are conducted alternatively (Park et al., 2014; Kim et al., 2016; Kim et al., 2018; Yang et al., 2020). Optical remote sensing also has the disadvantage that it cannot detect AOD for the regions covered by clouds or the regions with extremely high concentrations (Li et al., 2005; Nichol et al., 2010; Zhao et al., 2019). Such missing values can increase the data uncertainty in the analyses of the Earth’s environment (Youn et al., 2020). If necessary, a stable, gap-filled AOD image database can be constructed with the help of statistical techniques. To date, various studies have been carried out for the gap-filling of satellite AOD products using spatial statistical methods such as Kriging (Yu et al., 2011) and Bayesian approaches such as BMA (Bayesian Model Ensemble) (Singh et al., 2017) and BME (Bayesian Maximum Entropy) (Tang et al., 2016). AI (Artificial Intelligence) models were also built for gap filling of AOD using multiple explanatory variables related to AOD (Zhao et al., 2019; She et al., 2020). However, the approaches based on Bayesian statistics and AI models could not ensure a semi-real-time operational application for the gap-filling of hourly AOD products if an amount of auxiliary data for explanatory variables are not timely prepared.

Under the assumption of the preparation for realtime operation, this paper presents a spatial gap-filling framework using a univariate statistical method such as DCT-PLS (Discrete Cosine Transform-based Penalized Least Square Regression) (Garcia, 2010; Wang, 2012) and FMM (Fast Matching Method) inpainting (Telea, 2004). We conducted a feasibility test for the hourly AOD product from AHI (Advanced Himawari Imager) between January 1 and December 31, 2019, and compared the accuracy statistics of the two spatial gap-filling methods.

2. Data and methods

1) Himawari-8 AOD product

JMA (Japan Meteorological Agency) has launched Himawari-8, a next-generation geostationary meteorological satellite, on October 7, 2014, and are providing meteorological products such as AMV (Atmospheric Motion Vector), CSR (Clear Sky Radiance), HCAI (High-resolution Cloud Analysis Information), AOT, and ASWind (AMV-based Sea Surface Wind). The AHI onboard Himawari-8 has 16 spectral bands, including visible and infrared radiation (Table 1). The AOT or AOD product is created by the JAXA (Japan Aerospace Exploration Agency) algorithm (Yoshida et al., 2018) and provided in the format of NetCDF (Network Common Data Form) for the information of 500 nm AOD, AE (Ångström Exponent), and QA (Quality Analysis) flag. The spatial resolution is 0.05°, and the temporal resolution is 10 minutes for level 2 data; one hour, one day, and one month for level 3 data.

Table 1. Spectral information of AHI (Advanced Himawari Imager) (JMA, 2021)

The hourly product has two versions: AOT_Pure and AOT_Merged. We used the AOT_Merged because it is evaluated as a higher-quality product compared to AOT_Pure due to a statistical treatment considering the spatial and temporal variations of AOD (Kikuchi et al., 2018). This product is provided during the daytime when solar light is available (Lee and Lee, 2018). Hourly product has ten timeslots between 23 UTC (08 KST) and 09 UTC (18 KST) according to the daytime in Korea. We excluded two timeslots at 23 UTC close to dawn and 09 UTC close to sunset and used the eight images per day between 00 UTC and 07 UTC for our experiment.

2) DCT-PLS (Discrete Cosine Transform-based Penalized Least Square Regression)

DCT-PLS is a method devised for smoothing multidimensional data (Garcia, 2010; Wang, 2012). PLS (Penalized Least Square) regression pursuits a balance between original and smoothed data by minimizing Formula (1) that consists of a residual term between original and smoothed data and a penalty term for the roughness of the smoothed data (Whittaker, 1923; Wahba, 1990; Eilers, 2003).

\(F(\hat{y})=R S S+s P(\hat{y})=\|\hat{y}-y\|^{2}+s P(\hat{y})\)       (1)

where y is an original data; ˆyis the smoothed data; RSS is the residual sum of squares; sP(ˆy)is the penalty for the smoothed data. Because PLS can be formulated by the DCT (Discrete Cosine Transform) for multidimensional data (Garcia, 2010; Wang, 2012), DCT-PLS can be applied spatially for the gap-filling of gridded data.

\(F(\hat{X})=\left\|W^{1 / 2}{ }^{\circ}(\hat{X}-X)\right\|^{2}+s\|\Delta \hat{X}\|^{2}\)      (2)

where X is an original data with missing pixels; Xˆis the gap-filled data; Wis a weight matrix having the same dimension as X(a binary matrix with 0 for null pixels and 1 for not null pixels); | · | is Euclidean norm; Δ is Laplace operator; ° is Hadamard product; sis a smoothing parameter to overcome over or under-smoothing.

3) FMM (Fast Marching Method)

FMM inpainting technique, a weighted average using neighboring pixels around a null pixel, is appropriate for high-resolution image processing because of its fast calculation (Telea, 2004). If Ω is the null-pixel area to be filled and ε is the search radius, a temporal procedure is necessary to determine the sequence for Ω to fill from border to inside (Fig. 1). To fill the pixel p on the border of Ω, all the values q inside B_ε(p) are summarized by a normalized weighting function w(p, q) consisting of three components, namely, directional, geometric, and level-set weights. An image pixel value I(p) is defined as

Fig. 1. Principle of FMM inpainting (Telea, 2004).

\(I(p)=\frac{\sum_{q \in B_{ε}(p)} w(p, q)[I(q)+\nabla I(q)(p-q)]}{\left.\sum_{q \in B_{ε}(p)}\right) w(p, q)}\)       (3)

\(w(p, q)=\operatorname{dir}(p, q) * d s t(p, q) * \operatorname{lev}(p, q)\)       (4)

\(\operatorname{dir}(p, q)=\frac{p-q}{\|p-q\|} \cdot \mathrm{N}(\mathrm{p})\)       (5)

\(d s t(p, q)=\frac{d_{0}^{2}}{\|p-q\|^{2}}\)       (6)

\(\operatorname{lev}(p, q)=\frac{T_{0}}{1+|T(p)-T(q)|}\)       (7)

where I is the image value for a pixel q; Tis the propagated value for a pixel; ∇I denotes the image gradient; N denotes the normal direction that the propagation proceeds. dir(p, q) is the weight for the directional component; dst(p, q) is the weight for geometric distance; lev(p, q) is the weight for the level set distance that pixels close to the contour through p contribute more than farther pixels (Telea, 2004). Both dst(p, q) and lev(p, q) are relative with respect to the reference distance d0and T0, which indeed are set to the interpixel distance, i.e., to 1.

4) Experiment setup

The study area for a gap-filling experiment of AHI hourly AOD product is around South Korea (33.73~38.72°N and 125.78~129.77°E) with 8, 000 pixels (100×80) in a 0.05° grid. The experiment period is between January 1 and December 31, 2019. Out of the 2, 920 (365×8) images for daytime (00 to 07 UTC), we used 2, 914 except for the six missing images. Table 2 shows the null pixel statistics for the 2, 914 images. Based on the statistics, we applied DCT-PLS and FMM to only 201 images with a null pixel ratio under 0.6. For an objective performance test, we created two random blocks with 10×10 null pixels in the AOD images and conducted the spatial gap-filling followed by the pixel-to-pixel comparisons between the original and the gap filled blocks. If a null pixel takes up more than half of a random block, it is discarded, and a new random block is created. This produced 30, 164 pixels for the accuracy validation. Under such a setting, we carried out the spatial gap-filling of the AOD images using DCT-PLS and FMM and calculated the accuracy statistics for the 30, 164 pixels in terms of MBE (Mean Bias Error), MAE (Mean Absolute Error), RMSE (Root Mean Square Error), and CC (Correlation Coefficient).

Table 2. Null pixel statistics for AHI hourly AOD data, January 1 to December 31, 2019

3. Results and discussions

Regarding the 30, 164 pixels from the 201 images of AHI hourly AOD, we calculated the accuracy statistics by comparing the original and gap-filled images (Tables 3 and 4). DCT-PLS produced the CC of 0.988 and the MAE of 0.020. In the case of FMM with the search radius of 3, 4, 5, and 6 pixels, the result from ε=3 or ε=4 was slightly better than others. Because a too narrow or wide radius can increase the uncertainty in the result (Fan et al., 2013), ε=4 was determined as optimal in this experiment. The result of FMM had a CC of 0.980 and MAE of 0.028, which is very similar to that of DCT-PLS. When a gap-filling for a wider area is required, a wider search radius is necessary unless it worsens accuracy. Indeed, the setup of the search radius should be determined by the characteristics of the missing pixel distribution.

Table 3. Accuracy statistics of gap-filled AHI hourly AOD using DCT-PLS

Table 4. Accuracy statistics of gap-filled AHI hourly AOD using FMM

Fig. 2 is the scatterplot to show the comparison between the actual and gap-filled AOD values by DCT-PLS, and Fig. 3 is the gap-filling result of FMM (ε=4). The random null blocks for validation were filled with very similar values to the actual AOD by both DCT-PLS and FMM, so the scatterplots appeared concentrated around the 1:1 line.

Fig. 2. Observed vs. predicted AHI hourly AOD using DCT-PLS.

Fig. 3. Observed vs. predicted AHI hourly AOD using FMM.

Table 5 is the AOD statistics gathered from all pixels of the 201 images: (a) original images with random gaps, (b) gap-filled images using DCT-PLS, and (c) gap-filled images using FMM. Minimum values for the original image, DCT-PLS, and FMM were 0.00, and the maximum values were 2.37 for the three datasets. The mean value was 0.33 for the original data and 0.32 for DCT-PLS and FMM; the standard deviation was 0.13 for the original data and DCT-PLS and 0.14 for FMM. Also, the histogram for the three datasets was very similar, but a slight underestimation by PCT-PLS and FMM was shown (Fig. 4).

Table 5. Hourly AOD data statistics before and after gap-filling

Fig. 4. Hourly AOD histograms: (a) original images with random gaps, (b) gap-filled images using DCT-PLS, and (c) gap-filled images using FMM.

Fig. 5 shows a few examples for the map of gap filled AHI hourly AOD using DCT-PLS and FMM. In the case of March 19 (07 UTC), the null-pixel area is not large, and the two methods yield a similar result. The original data of June 4 (07 UTC) had a large null pixel area around the DMZ (demilitarized zone). The result of DCT-PLS showed a gradual change, but FMM produced an abrupt change in the AOD values because it fills the border pixel first and moves to the nearest neighbor. Overall, both methods showed fast calculation and high accuracy for the images with the null pixel ratio under 0.6. Also, they are the univariate gap-filling method that does not require auxiliary explanatory variables, which means that they can be applied to a real-time operational framework for the spatial gap-filling of satellite images.

Fig. 5. Maps for gap-filled AHI hourly AOD by DCT-PLS and FMM.

4. Conclusions

This paper examined the spatial gap-filling methods for the AHI hourly AOD product using DCT-PLS and FMM inpainting to solve missing values and conducted the feasibility tests for quantitative validation of the methods. When the null-pixel area is not very large (null pixel ratio < 0.6), both methods produced a high accuracy of the CC > 0.98 for the random blind tests. Since they are a univariate method with no need for additional explanatory variables, they could be used as a real-time operational algorithm for spatial gap-filling. However, most of the AHI hourly AOD images have a null pixel ratio over 0.6, so the other gap-filling methods using AI with extra explanatory variables may be necessary. These univariate and multivariate gap-filling methods should be tested for the low-resolution images from the AMI (Advanced Meteorological Imager) of GK2A (Geostationary Korea Multi-purpose Satellite 2A), GEMS (Geostationary Environment Monitoring Spectrometer) and GOCI2 (Geostationary Ocean Color Imager) of GK2B (Geostationary Korea Multi-purpose Satellite 2B) and the high-resolution images from the CAS500 (Compact Advanced Satellite) series soon.