Cloud-Type Classification for Southeast China Based on Geostationary Orbit EO Datasets and the LighGBM Model

Abstract

The study of clouds and their characteristics provides important information for understanding climate change and its impacts as it provides information on weather conditions and forecasting. In this study, Earth observation (EO) data from the FY4A AGRI and Himawari-8 CLP products were used to classify and identify distinct cloud types in southeastern China. To reduce the impact of parallax between geostationary satellites, we proposed adopting a sliding detection method for quality control of cloud-type data. Additionally, the Bayesian optimization method was employed herein to tune the hyperparameters of the LightGBM model. Our study results demonstrated that Bayesian optimization significantly increased model performance, resulting in successful cloud-type classification and identification. The simultaneous use of visible and shortwave infrared channels, and brightness temperature difference channels, enhanced the model’s classification performance. Those channels accounted for 43.79% and 21.84% of the overall features, respectively. Certainly, the model in this study outperformed compared with the traditional thresholding method (TT), support vector machine (SVM), and random forest (RF). Results showed a model prediction accuracy of 97.54%, which was higher than that of TT (51.06%), SVM (96.47%), and RF (97.49%). Additionally, the Kappa coefficient of the model was 0.951, indicating the model’s classification results were consistent with the true values. Notably, this performance also surpassed TT (0.351), SVM (0.929), and RF (0.950).

1. Introduction

Clouds have a significant role in the evolution of weather, climate, the global water cycle, and the global radiative budget with an average annual coverage of more than 60% of the Earth [1,2,3,4,5]. Different cloud types provide crucial cues for weather analysis and forecasting by reflecting various atmospheric states and circulation scenarios. Cloud types have clear indications of weather conditions and weather changes. For example, stratocumulus and nimbostratus are low clouds with clear indications of stratiform precipitation, cumulus, and cumulonimbus are effective indicators of convective precipitation, and cirrus clouds are generally not accompanied by precipitation [6]. In addition, the cloud classification products are important for the accuracy improvement of the retrieval results of meteorological parameters, such as cloud base height and precipitation intensity [7,8,9].

The southeastern region of China is located in the obvious monsoon zone. The abundant water vapor, as well as the power and thermal effects of the Tibetan Plateau, provide favorable conditions for the formation of clouds and precipitation [10,11,12]. The strength of the East Asian summer monsoon circulation affects the weather and climate of the region in significant ways. The seasonal rain belts, the beginning and end of the rainy season, and the structure of rain patterns are closely related to the East Asian summer monsoon [13]. Droughts and floods triggered by changes in regional precipitation seriously affect the production, life, and economic activities of local people. Furthermore, strong convective clouds may also trigger strong convective weather such as thunderstorms and gales, short-term heavy precipitation, and hail [14], resulting in serious economic losses and even casualties. Therefore, it is of great significance to achieve accurate cloud-type classification in the southeastern region of China.

Clouds are classified by meteorological operations into three families and ten genera based on their macro-structural features, including cloud base height, morphology, and structure. This method is considered the most effective way to distinguish cloud types artificially [6]. However, artificial identification has a certain degree of subjectivity, which will have a certain impact on the identification accuracy, and it is not possible to carry out continuous observation over a wide area. Geostationary satellites have the capability to carry out uniform and continuous observation of the same target area for a long time, making them an effective means of cloud observation. However, the traditional classification standard for ground-based cloud observations is not suitable for satellite data due to the different observation angles. The International Satellite Cloud Climatology Project (ISCCP) proposes a set of cloud classification standards [15], which classifies clouds into nine categories according to the cloud top pressure and optical thickness. The cloud classification product data specifically from the Himawari-8 geostationary satellite is based on the ISCCP standard.

Currently used cloud detection and classification algorithms can be classified into three categories: simple methods, statistical methods, and artificial intelligence methods [16]. The threshold-based classification is the most commonly used simple cloud classification method [17,18,19], which occupies less computational resources but has low accuracy. Statistical methods have higher accuracy compared to simple methods. However, those methods also require more computational power, such as K-means [20] and maximum-likelihood estimation [21]. With the recent technological developments, more and more research is focusing on the use of artificial intelligence in mapping clouds and their properties. Artificial intelligence is a method of using computers and machines to mimic the way the human brain thinks to solve problems and make decisions, and machine learning is a subfield of artificial intelligence. Traditional machine learning requires manual extraction of feature parameters for model training. Classical machine learning-based methods train models to classify cloud types by manually selecting feature parameters, such as brightness temperature, and texture. Such methods, for example, Support Vector Machines [22] and random forests [21,23], have demonstrated that machine learning methods can be successfully used for cloud detection and classification tasks.

The FY-4A geostationary satellite has no product based on the ISCCP cloud classification standard, but some researchers have already conducted research on FY-4A based on the cloud-type classification product of Himawari-8 [24,25]. Nevertheless, they mostly mapped the Himawari-8 satellite data directly to FY-4A data without considering the inter-satellite parallax effect, which may lead to the wrong matching of cloud-type labels. For supervised learning, label mislabeling will have a great impact on the model accuracy. Therefore, it is an imperative step to perform quality control on labeled data before constructing the dataset to minimize matching errors and the impact of data quality on the model.

Among machine learning algorithms, LightGBM has become a research hotspot in recent years due to its fast training speed, low resource consumption, strong generalization ability, and it is not prone to overfitting. LightGBM is a decision tree-based gradient-boosting algorithmic framework proposed by Microsoft in 2017 [26]. It is widely used in various fields such as geomatics [27], Earth Observation (EO) [28], and meteorology [29].

In the purview of the above, the present study aimed to propose a cloud-type classification retrieval based on the LightGBM model. Currently, there have been several studies on the development of AGRI cloud classification products, but most of them do not consider the effect of parallax. This is where the contribution of this study lies. A sliding detection method was used for the quality control of the cloud-type data and to reduce the effect of parallax on the model. Bayesian optimization was used to optimize the model hyperparameters. Furthermore, comparisons were made with other cloud classification methods to assess the accuracy of the retrieval model in this paper.

2. Data and Methods

2.1. Data

2.1.1. FY-4A Dataset

FY-4A is a Chinese new-generation geostationary meteorological satellite, launched on 11 December 2016. The on-board Advanced Geostationary Radiation Imager (AGRI) includes 14 channels capable of providing minute-level observations. Compared with the previous generation of geostationary meteorological satellites, the Feng Yun 2 (FY-2) series of VISSR, AGRI has greatly improved in terms of spatial and temporal resolution and the number of observation channels [30]. The observation parameters of each channel are shown in

Table 1. FY-4A AGRI channel information.

Band	Spectral Coverage	Central Wavelength	Spectral Bandwidth	Main Applications
1	VIS/NIR	0.47 µm	0.45~0.49 µm	Aerosol, visibility
2		0.65 µm	0.55~0.75 µm	Fog and cloud
3		0.825 µm	0.75~0.90 µm	Vegetation
4	Shortwave IR	1.375 µm	1.36~1.39 µm	Cirrus cloud
5		1.61 µm	1.58~1.64 µm	Cloud and snow
6		2.25 µm	2.1~2.35 µm	Cirrus clouds and aerosol
7	Midwave IR	3.75 µm	3.5~4.0 µm	Fire point
8		3.75 µm	3.5~4.0 µm	Earth’s surface
9	Water vapor	6.25 µm	5.8~6.7 µm	Upper-level WV
10		7.1 µm	6.9~7.3 µm	Mid-level WV
11	Longwave IR	8.5 µm	8.0~9.0 µm	Cloud motion wind and cloud
12		10.7 µm	10.3~11.3 µm	Sea surface temperature
13		12.0 µm	11.5~12.5 µm	Sea surface temperature
14		13.5 µm	13.2~13.8 µm	Cloud top height

.

The AGRI full disk data from June to August 2022 were used in this paper. The temporal resolution of the data was 15 min and the spatial resolution was 4 km. The data were provided by the National Meteorological Satellite Centre (http://satellite.nsmc.org.cn/ (accessed on 13 March 2023)). Since the Himawari-8 CLP product only has data during the daytime hours, we only selected the data for the whole point and half point, with a time range of 00-09UTC. We set the southeast coastal area of China as the study area, with a longitude and latitude range of 105°~123°E, 15°~35°N. The study area is shown in Figure 1.

2.1.2. Himawari-8 Dataset

The Himawari-8 satellite, launched on 7 October 2014, is a new-generation geosynchronous satellite developed by the Japan Aerospace Exploration Agency (JAXA). The Advanced Himawari Imager (AHI) is one of the main payloads of the Himawari-8 satellite, with a spatial resolution of 0.5–2 km and a temporal resolution of 10 min. The AHI sensor has 16 bands and provides observational support for Japan, East Asia, and the entire western Pacific region with a temporal resolution of 10 min. JAXA generates a wide range of geophysical observations using Himawai-8 standard data, and the cloud-type data used in this paper were provided by JAXA’s P-tree system (https://www.eorc.jaxa.jp/ptree/index.html (accessed on 23 March 2023)) for model parameter acquisition and accuracy validation. The CLP product consisted of parameters such as cloud optical thickness, cloud effective particle radius, cloud top temperature, cloud top height, cloud type defined by ISCCP, and quality control information, with a data resolution of 0.05°, which is in good agreement with the MODIS product [31]. According to the ISCCP standards, the cloud types included cirrus (Cirrus, Ci), cirrostratus (Cs), deep convection (Deep Convection, Dc), altocumulus (Altocumulus, Ac), altostratus (Altostratus, As), and nepheloid (Nimbostratus, Ns), Cumulus (Cu), Stratocumulus (Sc), and Stratus (St) in nine categories, as shown in Figure 2. The Himawari-8 satellite cloud-type product is shown in Figure 3.

2.2. Data Pre-Processing

2.2.1. AGRI Data Preprocessing

In order to facilitate the calculation, data pre-processing of the L1 data was required. Firstly, the data was radiometrically calibrated, and DN values were converted to either albedo or cloud top brightness temperature based on the CAL Channel calibration table. The FY-4A satellite uses a nominal projection of the geostationary orbit as defined by the CGMS LRIT/HRIT global specification, which requires the calculation of latitude and longitude based on the WGS84 reference ellipsoid after radiometric calibration. After that, the data were interpolated onto a 0.05° × 0.05° grid of equal latitude and longitude, which was consistent with the resolution of the CLP product. Different light source conditions at different moments can cause inconsistencies in the brightness of the visible (VIS) image, and to eliminate this effect, the sun altitude angle was revised to ensure that the radiance of the pixels is consistent at different sun altitude angles. In addition, the data ranges of each channel were different, and to facilitate the training process, the data of different channels needed to be normalized by using the maximum-minimum processing method so that all the data were between 0 and 1. The calculation method i as follows:

$${\mathrm{x}}^{\prime }=\frac{\mathrm{x}-\min }{\max -\min }$$

(1)

where ${\mathrm{x}}^{\prime }$ is the data obtained after normalization, x is the original data of each channel, min is the minimum value of each channel, and max is the maximum value of each channel.

2.2.2. Cloud Type Label Matching

The parallax effect needed to be taken into account when observations from instruments on both geostationary satellites, FY-4A and Himawari-8, are used jointly. Parallax is the deviation between the satellite observation target position and the true position, it is prevalent in non-vertical observations [33]. When the satellite observation target has a large viewing angle, the cloud target position measured by the satellite is not the true cloud position, but will be shifted in the direction away from the sub-satellite point, and this error is especially noticeable for samples far away from the sub-satellite point. As shown in Figure 4, When both AGRI and AHI observe the same target position, point P, O1 represents the observation cloud position for AGRI, while O2 represents the observation cloud position for AHI. The parallax between the two satellites is represented by d. In this case, the radiation obtained from the observation at point P corresponds to the observation at point O1, and the cloud-type label obtained from the observation at point P corresponds to the observation at point O2. To reduce the possible matching errors, the CLP product needed to be quality controlled by considering only the samples inside the cloud body and ignoring the samples at the edges to reject the incorrectly matched cloud labeling data.

In Figure 4, α and θ are the zenith and azimuth angles of AGRI, β, and γ are the zenith and azimuth angles of AHI, a is the parallax of AGRI, b is the parallax of AHI, d is the parallax between the two satellites, and h is the cloud height. Based on the latitude and longitude of O1 and O2 and the position of the satellite’s subsatellite point, the zenith angle and azimuth of each point can be calculated. With the known cloud height, d can be calculated [33]:

$$\mathrm{a}=\mathrm{h}\ast \tan \mathsf{\alpha }$$

(2)

$$\mathrm{b}=\mathrm{h}\ast \tan \mathsf{\beta }$$

(3)

$$\angle \mathrm{O}1{\mathrm{P}}^{\prime }\mathrm{O}2=\mathsf{\gamma }-\mathsf{\theta }$$

(4)

$$\mathrm{d}=\sqrt{{\mathrm{h}}^2\left({\tan }^2\mathsf{\alpha }+{\tan }^2\mathsf{\beta }\right)-2{\mathrm{h}}^2\tan \mathsf{\alpha }\tan \mathsf{\beta }\cos (\mathsf{\gamma }-\mathsf{\theta })}$$

(5)

It can be observed from the equation that inter-satellite parallax depends on the zenith angle, azimuth angle, and cloud top height. When FY-4A and Himawari-8 observe jointly, the parallax between the two satellites can be calculated by knowing the cloud top height, satellite zenith angle, and satellite azimuth angle. The cloud top heights in the southeast region of China are relatively low, mostly medium-low clouds, and the average cloud top heights are mostly between 7–10 km [30]. According to Equations (2)–(5), the actual parallax between the two satellites can be up to 12 km (about 3 pixels), and thus the 7 × 7-pixel matrix sliding detection method was used to cull the data. The step size was 1. When all cloud types within 7 × 7 pixels were consistent, it was considered that the cloud type label of the center pixel could be matched correctly; otherwise, the data of this center pixel point was rejected.

Ensuring the temporal and spatial correspondence between the input feature variables and the real target is necessary during model training and testing. Therefore, after quality control of cloud-type data, we need to perform spatio-temporal matching between AGRI images and CLP products. Through data preprocessing, we reprojected the AGRI data onto a 0.05 × 0.05° grid of equal latitude and longitude, which is consistent with the spatial resolution of the CLP product. Therefore, we selected data files with the same observation time and directly mapped the processed CLP product onto the AGRI image to complete the data matching. Taking 1 July 2022 at 06:00 UTC as an example, Figure 5 shows the cloud types and corresponding AGRI cloud images for the ISCCP standard after time-space matching. Additionally, we need to remove pairs with incomplete features and missing cloud-type labels to construct the dataset.

Ultimately, the AGRI and AHI dataset covering June–August 2022 was created, labeling over 9,900,000 pixels. We randomly selected 80% of the data in the dataset as the training set and the remaining 20% as the testing set. They were used for model training as well as validation, respectively.

2.3. LightGBM

LightGBM is a distributed boosting algorithm framework based on the decision tree algorithm, which is an efficient implementation of the gradient boosting decision tree (GBDT) algorithm. GBDT is widely used for classification, regression, sorting, and other tasks due to its high efficiency, accuracy, interpretability, and strong generalization ability. It is applicable to various classification problems and is not prone to overfitting. It takes the negative gradient of the loss function as a residual approximation of the current decision tree, goes on to fit a new decision tree, and accumulates the results of multiple decision trees as the final predicted output. The objective function can be expressed as:

$${\mathrm{Obj}}^{\left(\mathrm{t}\right)}=\sum _{\mathrm{i}=1}^{\mathrm{n}}\mathrm{L}({\mathrm{y}}_{\mathrm{i}},{\widehat{{\mathrm{y}}_{\mathrm{i}}}}^{(\mathrm{t}-1)}+{\mathrm{f}}_{\mathrm{t}}({\mathrm{x}}_{\mathrm{i}}\left)\right)+\sum _{\mathrm{i}=1}^{\mathrm{t}}\mathsf{\Omega }({\mathrm{f}}_{\mathrm{i}})$$

(6)

$$\mathsf{\Omega }({\mathrm{f}}_{\mathrm{i}})=\mathsf{\gamma }\mathrm{T}+\mathsf{\lambda }\frac{1}{2}\sum _{\mathrm{j}=1}^{\mathrm{T}}{\mathrm{w}}_{\mathrm{j}}^2$$

(7)

where L is the loss function, which measures the difference between the true value $y_i$ and the predicted value, ${\widehat{\mathrm{y}}}^{(\mathrm{t}-1)}$ is the result of the (t − 1)-th iteration, $\mathsf{\Omega }({\mathrm{f}}_{\mathrm{i}})$ is a regular term, γ can control the number of leaf nodes, T is the number of leaf nodes, and w is the number of leaf nodes. The objective function is to find a suitable tree ${\mathrm{f}}_{\mathrm{t}}$ that minimizes the value of the function.

We used the Taylor formula to expand the objective function to the second order:

$${\mathrm{Obj}}^{\left(\mathrm{t}\right)}=\sum _{\mathrm{i}=1}^{\mathrm{n}}\left[\mathrm{L}\right({\mathrm{y}}_{\mathrm{i}},{\widehat{{\mathrm{y}}_{\mathrm{i}}}}^{(\mathrm{t}-1)})+{\mathrm{g}}_{\mathrm{i}}{\mathrm{f}}_{\mathrm{t}}({\mathrm{x}}_{\mathrm{i}})+\frac{1}{2}{\mathrm{h}}_{\mathrm{i}}{\mathrm{f}}_{\mathrm{t}}^2({\mathrm{x}}_{\mathrm{i}}\left)\right]+\sum _{\mathrm{i}=1}^{\mathrm{t}}\mathsf{\Omega }({\mathrm{f}}_{\mathrm{i}})$$

(8)

$${\mathrm{g}}_{\mathrm{i}}=\frac{\partial \mathrm{L}({\mathrm{y}}_{\mathrm{i}},{\widehat{{\mathrm{y}}_{\mathrm{i}}}}^{(\mathrm{t}-1)})}{\partial {\widehat{{\mathrm{y}}_{\mathrm{i}}}}^{(\mathrm{t}-1)}}$$

(9)

$${\mathrm{h}}_{\mathrm{i}}=\frac{{\partial }^2\mathrm{L}({\mathrm{y}}_{\mathrm{i}},{\widehat{{\mathrm{y}}_{\mathrm{i}}}}^{(\mathrm{t}-1)})}{\partial {\widehat{{\mathrm{y}}_{\mathrm{i}}}}^{(\mathrm{t}-1)}}$$

(10)

With the problem that traditional GBDT causes a huge overhead in time and space, which is caused by the pre-sorting algorithm to find the data splitting points, LightGBM proposes the gradient-based one-sided sampling (GOSS) and exclusive feature bundling (EFB) techniques. The techniques optimize the data with large sample sizes and high feature dimensions, respectively, and prove on several public datasets that LightGBM improves the training speed by more than 20 times while achieving similar accuracy to GBDT [26]. LightGBM uses a Leaf-wise growth strategy to find one leaf at a time with the largest splitting gain from all the current leaves and then split it. With the same number of splits, Leaf-wise can reduce more error than Level-wise and have better accuracy. The schematic is shown in Figure 6. In this paper, we proposed a cloud classification model based on the LightGBM algorithm, implemented with Python using lightgbm package version 3.3.5.

2.4. Bayesian-Optimization

Bayesian optimization is a powerful strategy for finding the extremes of an objective function [34,35]. When searching for parameters that can maximize the objective function globally, existing prior information is taken into account to better tune the current parameters. Compared with other hyperparameter optimization schemes, like grid search and random search, Bayesian optimization is faster and can effectively improve the efficiency of finding hyperparameters. In this paper, we optimized the parameters of the LightGBM algorithm under the Bayesian framework.

Bayesian optimization uses a Gaussian process to fit the optimization objective function with reference to the previous parameter information and continuously updating the prior knowledge. The objective function for the hyperparametric optimization of the LightGBM model can be expressed as [36]:

$${\mathrm{x}}^{\ast }=\underset{\mathrm{x}\in \mathsf{\chi }}{\mathrm{argmax}}\mathrm{f}\left(\mathrm{x}\right)$$

(11)

where x is the hyperparameters of the model to be optimized, $\mathsf{\chi }$ is the search space of the hyperparameters, and f(x) is the objective function between the hyperparameters and the model performance. Using Bayesian optimization aims to find the appropriate combination of hyperparameters to maximize the model performance. The general procedure of Bayesian optimization for iterative search of hyperparameter combinations is as follows:

The loop iterates t times (t = 1,2,3…):

• According to the maximized acquisition function $\mathsf{\alpha }\left(\mathrm{x}\right)$, to find the next set of possible evaluation points, ${\mathrm{x}}_{\mathrm{t}+1}=\underset{\mathrm{x}\in \mathsf{\chi }}{\mathrm{argmax}}\mathsf{\alpha }(\mathrm{x};{\mathrm{D}}_{1:\mathrm{t}})$.
• Calculate the corresponding model performance based on the evaluation points ${\mathrm{x}}_{\mathrm{t}+1}$.
• Update ${\mathrm{x}}_{\mathrm{t}+1}$ to the previous observation ${\mathrm{D}}_{1:\mathrm{t}+1}=\{{\mathrm{D}}_{1:\mathrm{t}},({\mathrm{x}}_{\mathrm{n}+1},{\mathrm{f}}_{\mathrm{n}+1}\left)\right\}$, and update the probabilistic surrogate model.

The loop ends and the final combination of hyperparameters.

2.5. Evaluation Indicators

Several previous studies have used the same metrics to evaluate the classification model, including Precision (P), Recall (R), F1 Score (F1), Accuracy (Acc), and Kappa Coefficient (Kappa) [37,38,39,40]. The formula for each metric is as follows:

$$\mathrm{P}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}$$

(12)

$$\mathrm{R}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$$

(13)

$$\mathrm{F}1=2\times \frac{\mathrm{P}\times \mathrm{R}}{\mathrm{P}+\mathrm{R}}$$

(14)

$$\mathrm{Acc}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}$$

(15)

$$\mathrm{Kappa}=\frac{\mathrm{N}\sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{x}}_{\mathrm{ii}}-\sum _{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{x}}_{\mathrm{i}+}{\mathrm{x}}_{+\mathrm{i}})}{{\mathrm{N}}^2-\sum _{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{x}}_{\mathrm{i}+}{\mathrm{x}}_{+\mathrm{i}})}$$

(16)

TP is the number of pixels whose true sample category is class a and the recognition result is also class a, TN is the number of pixels whose true sample category and recognition result are both other categories, FP is the number of pixels that should have been identified as class a but were incorrectly identified as other categories, and FN is the number of pixels that should have been identified as other categories but was incorrectly identified as class a. N is the total number of samples, ${\mathrm{x}}_{\mathrm{ii}}$ is the number of correctly recognized sample size, and ${\mathrm{x}}_{\mathrm{i}+}$, ${\mathrm{x}}_{+\mathrm{i}}$ are the sum of the sample size of each row and column of the confusion matrix, respectively.

The above precision, recall, and F1 scores are only applicable to binary classification problems. For multiclassification problems, the Macro Average can be used as an evaluation metric to assess the classification performance of the model in all classes. Macro Average is the arithmetic mean of the performance metrics for each category, which can clearly show the recognition effect of different categories of samples under the unbalanced state.

2.6. Experimental Setup

The flowchart in Figure 7 illustrates the main content of this paper. It is mainly divided into two parts: data preprocessing and model training and evaluation. The data preprocessing stage, including the preprocessing of AGRI data, the quality control of CLP products, and the dataset construction, will be described in detail in this section. The model construction and evaluation will be given in Section 3.

In order to find the optimal channel combination for cloud classification, we designed three sets of plans to verify the effects of different channel combinations on the experimental results, and the specific channel combination plans are shown in

Table 2. Different channel combination settings.

No.	Channel Combinations
1	Chn07~Chn14
2	Chn01~chn14
3	Chn01~chn14+ BTD (10−12) + BTD (11−12) + BTD (12−13)

. The input data of Plan 1 was 7–14 channel brightness temperature, the input data of Plan 2 was 1–14 channel brightness temperature, and Plan 3 added the BTD (10–12), BTD (11–12), and BTD (12–13) three brightness temperature difference, commonly used in cloud analysis [19,39]. BTD (10–12) was used to indicate the height of the cloud development, BTD (11–12) was used to differentiate between water and ice clouds, and BTD (12–13) was the split window channel brightness temperature difference that was used to distinguish between thin and thick clouds.

3. Results

3.1. Model Tuning

In the present study, Bayesian optimization was used to tune the hyperparameters of the model, and the area under the curve (AUC) was employed as the performance metric [41], which was insensitive to class skew. The optimizer used a 5-fold cross-validation scheme, and the optimal hyperparameter solution for the model was finally obtained by performing several iterations of the optimization process.

The hyperparameter optimization process corresponding to the three plans is shown in Figure 8, where the blue points indicate the results of each iteration, the orange points indicate the optimal results during the iteration process, and the labeled information is the corresponding optimal hyperparameter combination. From the figure, it can be observed that during the process of hyperparameter tuning, the AUC of each experimental group was improved to a certain extent. The optimum was reached in the 16th round of parameter tuning for Plan 1, which was 0.975, the highest value of cross-validated AUC for Plan 2 was 0.973, which was obtained in the 39th round of tuning, and, for Plan 3, the corresponding highest value was 0.994, which was reached in the 44th round of tuning.

We retrained the models using the optimal hyperparameter combinations from the three sets of plans and obtained three models: Model #1, Model #2, and Model #3. To show more intuitively the improvement of the training effect of hyperparameter optimization on the models, the training curves before and after hyperparameter optimization are compared in Figure 9. As can be seen from the figure, when the models were trained with the default hyperparameters, the log loss curves showed an obvious and drastic jitter, the accuracy curves of the models were also unstable, and the training curves of the three models showed different degrees of jitter after the 54th, 48th, and 43rd iterations, respectively. In contrast, the three models using the optimal hyperparameter combinations had good fitting results, with Acc on the training set reaching 0.921, 0.973, and 0.975, respectively.

3.2. Comparison of Different Channel Combination Models

In order to quantitatively assess the classification results of different channel combinations, the classification results of the three models were compared pixel-by-pixel with the quality-controlled Himawari-8 CLP product and evaluated according to different evaluation metrics. The classification results of the three models are shown in Figure 10. The comparison results of different channel combinations on the test set showed that the Acc of Model #2 improved from 92.10% to 97.33% as compared to Model #1, but was still lower than 97.54% of Model #3. Model #1 did not recognize Ac and Cu well, with a precision of 19.05% and 68.75%, while recall was only 2.31% and 0.91%. In contrast, Model #2 had a great improvement in recall for these two types of clouds, reaching 31.21% and 61.10%, respectively, and the F1 score was still higher than that of Model #1 despite a certain decrease in the precision rate, which was mainly due to the advantage of VIS channels in recognizing the thickness of the cloud. For different types of clouds, Model #3 always had the highest precision, recall, and F1 score. The Kappa coefficients of the three models were 0.843, 0.947, and 0.951, and there was a high degree of agreement between the model classification results and the true values. Overall, Model #3 had the best classification results, so Model 3 was selected as the final cloud classification model.

3.3. Model Evaluation

Comparing different combinations of input channels, we found that for most cloud types, Model #3 was the best model. However, different feature variables contributed differently to the cloud classification task, and the importance of different channels needed to be analyzed. Figure 11 shows the order of importance of the features of the model, where a larger value indicates that the feature was more important. The vertical coordinates are the different characterizing variables of the model, and the horizontal coordinates are the scores of the characterizing variables.

The feature importance analysis showed that the top three feature variables among all the feature variables of the model were the first channel (0.47 µm) albedo, the BTD (12–13), and the BTD (11–12), which indicated that there were more significant differences in the values of these features among different cloud types. The VIS and Shortwave infrared features and the bright temperature difference features accounted for about 43.79% and 21.84% of the overall features. The ISCCP cloud classification standard classifies clouds based on the cloud optical thickness and the cloud top pressure, so the first-channel albedo, which reflects the thickness of the clouds, and the BTD (12–13), which characterizes the cloud optical thickness, contributed the most to the classification of the cloud-type. The phases of clouds at different heights are also different, with low clouds dominated by water clouds, medium clouds mainly ice-water mixed clouds, and high clouds ice clouds [32,42]. Therefore, the BTD (11–12), which can characterize the cloud phase, was also important.

A visualization of the confusion matrix of the model’s recognition results on the testing set (Figure 12). As can be observed, the first column is relatively darker, which indicates that the model preferred to recognize other types of clouds, especially those with thin cloud thickness (Ci, Ac, and Cu) as clear. Ci, Ac, and Cu are usually located at the edge of the cloud mass and make up a lower proportion of the training samples than the other types of clouds, resulting in an imbalance between the multiple sample categories. At this point, the model was not able to learn enough from a small number of sample categories, and the model preferred clear skies when distinguishing similar spectral features. In addition, since edge clouds are more susceptible to parallax, the quality-controlled data may still have label-matching errors, affecting the model’s classification accuracy.

3.4. Comparison of Model Accuracy

To further evaluate the cloud classification results achieved by the algorithms in this paper, a comparison was made with the traditional thresholding method (TT) [18,19], support vector machine (SVM), and random forest (RF). Figure 13 shows the AGRI 13-channel bright temperature images and the cloud-type classification results for the four algorithms. All four algorithms could detect clouds with high cloud top heights, as shown in the figure. It was difficult to distinguish between clear and low clouds using the TT algorithm because it lacked information about the VIS channels.

To accurately assess the classification results of various algorithms, we performed a quantitative comparison of recognition results. The results of each algorithm on the testing set are shown in the table below.

Table 3. Comparison of classification results of different algorithms.

Algorithms	Macro P	Macro R	Macro F1	Acc	Kappa
LightGBM	84.10%	78.24%	79.74%	97.54%	0.951
TT	26.60%	31.20%	22.92%	51.06%	0.351
SVM	70.00%	56.98%	59.63%	96.47%	0.929
RF	73.11%	73.03%	72.97%	97.49%	0.950

shows that TT had the lowest recognition accuracy. The accuracy of all machine learning models, except for the TT model, was above 90%, and the kappa coefficient was higher than 0.8, indicating a strong consistency between the model classifications and actual values. Overall, these three machine learning models showed good classification performance on the testing set. The LightGBM model had the highest performance on all evaluation metrics, reaching 84.10%, 78.24%, 79.74%, 97.54%, and 0.951, respectively. Combining all the evaluation indicators, the LightGBM model had the best performance in terms of results.

4. Discussion

In this study, a machine learning method was proposed that aimed to achieve a cloud classification of FY-4A satellite images (ISCCP standard) by introducing Bayesian optimization. Performing quality control on cloud-type data effectively improved data quality and reduced the effect of inter-satellite parallax on model accuracy. In addition, the optimal model for cloud-type classification was found, through a set of sensitivity tests of input channel combinations. Compared with many other excellent cloud classification methods, the model proposed in this paper also achieved the best results. However, the model still has some limitations in recognizing certain types of clouds. Our work has some implications for future research on geostationary satellite cloud-type classification:

• By comparing the performance of the model before and after hyperparameter optimization for the three different input channel combinations, it was demonstrated that Bayesian optimization effectively improved the model performance.
• In evaluating the performance of the different input channel combination models on the testing set, we found that the advantage of the VIS and Shortwave channels in recognizing cloud thickness improved the classification effect of the models to a certain extent, especially the recall of Ac and Cu, which are two types of cloud types with thinner cloud thickness and lower cloud top height. However, because of the spectral similarity between clear and Ac and Cu, caused the model to misclassify clear as these two types of clouds, leading to decreased precision.
• In evaluating the optimal model, we found that the VIS and Shortwave infrared channels and the brightness temperature difference channel accounted for 43.79% and 21.84% of the overall features, and the top three feature variables in terms of feature importance were the first-channel (0.47 μm) albedo, BTD (12–13), and BTD (11–12). When the model classifies cloud types, it is easy to identify clouds with thinner cloud thickness (Ci, Ac, and Cu) as clear. One reason for this is that Ci, Ac, and Cu clouds tend to be located at the edges of the cloud, which are represented by a smaller number of training samples than other cloud types. This can lead to the model not having enough learning ability for these categories. Additionally, edge clouds are more likely to be impacted by parallax effects, and, even after quality control, there may still be labeling errors that affect the accuracy of the model’s classifications.

5. Conclusions

In this paper, we proposed a cloud classification model for FY-4A satellite observation images that can be used in the southeast region of China. To address the problem of parallax between geostationary satellites, we used the sliding detection method for the quality control of the cloud-type label data to reduce the impact of label-matching errors on the accuracy of the model. The setting of hyperparameters will largely affect the classification effect of the machine learning model, and by introducing the Bayesian optimization method, we obtained the optimal hyperparameter combination of the model. We found that Bayesian optimization adjusted the model hyperparameters well. Compared with the model using only the bright temperature of the infrared channels, the introduction of albedo and bright temperature difference improved the sensitivity of the model for cloud thickness identification to some extent. Finally, we achieved a quantitative assessment of regional cloud-type classification. For validation, we compared it with the Himawari-8 CLP product after quality control and the Acc of the optimal model reached 97.54%, which was greater than the 51.06% of TT, 96.47% of SVM, and 97.49% of RF, and the model achieved the best results. Combined with the feature importance analysis, it further revealed the degree of influence of different input channels on the cloud classification task in the region. The results of this study can provide useful information for meteorological operations in the region.

However, some cloud types were not well recognized due to factors such as the model feature variables and the sample imbalance of each type of cloud. We plan to solve these problems by further improving the dataset and optimizing the model in our future work.