A knowledge-based data-driven (KBDD) framework for all-day identification of cloud types using satellite remote sensing

The satellites Himawari-8/9 belong to the Himawari 3rd generation programme [Bessho et al., 2016, Kurihara et al., 2016, Yang et al., 2020], whose main mission is operational meteorology and additional mission is environmental applications. The Himawari-8/9 satellites were launched on October 7, 2014 and November 2, 2016, respectively. The satellite Himawari-8 replaced the satellite MTSAT-2 as an operational satellite on July 7, 2015, and the satellite Himawari-9 replaced the satellite Himawari-8 as the primary satellite on December 13, 2022. The scope of the satellites Himawari-8/9 is shown in Fig. 1, whose coverage area is from 80°E to 160°W, and from 60°N to 60°S. Satellite spectral channels provide albedo from B01 to B06 and brightness temperature (BT) from B07 to B16, whose specific information is shown in Table 1 [Huang et al., 2017, Taniguchi et al., 2022]. In this study, Himawari L1 gridded data within the study area were downloaded through the Japan Aerospace Exploration Agency (JAXA) Himawari Monitor P-Tree System (available at https://www.eorc.jaxa.jp/ptree/index.html).

The remote sensing data of the spatial resolution $0.05^{\circ}\times 0.05^{\circ}$ at 03:00 (UTC+0) were downloaded every 3 d from January 1, 2020 to December 31, 2022. The value range of albedo from B01 to B06 is [0, 100]%, but the value range of BT from B07 to B16 needs to be determined statistically. The BT values corresponding to different channels of each pixel in the downloaded remote sensing images are counted. The frequency and cumulative percentage curves of BT from B07 to B16 are plotted in Fig. 2(a) and 2(b), respectively. In Fig. 2(a), the BT value of the peak frequency of B08 is the smallest, while the BT value of the peak frequency of B07 is the largest. The frequency curves for B11, B13, B14, and B15 have the same trend; those for B08, B09, and B10 are relatively consistent with a normal distribution; and the BT values for B12 are evenly distributed between 235K and 275K. In order to accurately determine the range of BT values for each channel, the cumulative percentage of BT for each channel is calculated in Fig. 2(b). The images at the 0% and 100% positions of the cumulative percentage curve are magnified and observed to determine the range of BT values for each channel. The range of BT values for each channel is determined and recorded in Table 1, which covers 99.9% of its own data. Based on the range of spectral channel values, the model inputs can be reasonably normalized.

The International Satellite Cloud Climatology Project (ISCCP) definition of cloud types is presented in Fig. 3, which is detailed at https://isccp.giss.nasa.gov. Under cloudy skies, cloud types are classified as cirrus (Ci), cirrostratus (Cs), deep convection (Dc), altocumulus (Ac), altostratus (As), nimbostratus (Ns), cumulus (Cu), stratocumulus (Sc), and stratus (St), depending on cloud optical thickness and cloud top pressure. Otherwise, the sky is clear (Cl). The JAXA Himawari Monitor P-Tree System provides Himawari L2 gridded data cloud types with spatial resolution $0.05^{\circ}\times 0.05^{\circ}$ and temporal resolution 10 min [Letu et al., 2019, Lai et al., 2019, Letu et al., 2020].

In Fig. 4, the change process of the overall cloud-type distributions every hour on June 10, 2021 is shown. Cloud types are found to be available only in the daytime region, while those in the nighttime region are unknown. Wang et al. [2022b] pointed out that the physical algorithms [Nakajima & Nakajma, 1995, Ishida & Nakajima, 2009] used by JAXA’s cloud-type product involve visible (VIS) B01 $\sim$ B03, near-infrared (NIR) B04 $\sim$ B06 and infrared (IR) B07 $\sim$ B16 bands, and thus cloud-type retrieval is limited to the daytime region only. Here, the cloud-type distributions at 03:00 (UTC+0) every 3 d are chosen as classification labels for this research from January 1, 2020 to December 31, 2022. Therefore, Himawari L1 gridded data with the same resolution at the corresponding time are downloaded in section 2.1.

The frequency of each cloud type is calculated, as shown in Fig. 5. The frequency of clear sky (Cl) in the entire satellite observation area is the highest, while the frequency of St is the lowest. Cl accounts for 26.42%, St accounts for 1.73%, and their ratio is approximately 15:1. The proportion of all types of clouds is 73.58%.

In order to derive all-day identification of cloud types through spectral information from Himawari-8/9 satellite sensors, a knowledge-based data-driven (KBDD) framework is proposed, whose specific architecture is depicted in Fig. 6. The KBDD framework mainly consists of knowledge module, mask module, addition of auxiliary information, network candidate set, and mask loss.

The main function of the knowledge module is to reorganize spectral information from Himawari-8/9 satellite sensors based on existing research knowledge. From previous research knowledge, it has been found that cloud types are not only related to single channel remote sensing, but also to the differences between two different channels. Here, it should be noted that satellite spectral channels B01 to B06 are characterized by albedo, while satellite spectral channels B07 to B16 are characterized by brightness temperature. Since albedo and brightness temperature are distinct variables and their units are different, the pairwise combinations must have the same dimension in order to make a difference. The knowledge module completed the reorganization of the satellite spectral information of B01 $\sim$ B16, the differences between the pairwise combinations of B01 $\sim$ B06, and the differences between the pairwise combinations of B07 $\sim$ B16.

The mask module is a convenient way to mask certain spectral channels and increase the flexibility of the framework, which has the ability to train with/without VIS and NIR data but does not require changing the model structure. When not using VIS and NIR channels, the mask module can set the satellite spectral information of B01 $\sim$ B06 and the differences between the pairwise combinations of B01 $\sim$ B06 to zero through the parameters mask ratio and mask bands.

The addition of auxiliary information is to improve the performance of the framework. The auxiliary information includes satellite zenith angle (SAZ), satellite azimuth angle (SAA), solar zenith angle (SOZ), solar azimuth angle (SOA), longitude, latitude, altitude, and underlying surface attributes (land or water). The auxiliary information and the reorganized satellite spectral information are merged through concatenation. The merged data is input into the network candidate set to obtain the probability result of cloud types. Due to the fact that the cloud-type references in the nighttime area do not exist, the mask loss mainly calculates the loss of labeled pixels in the target image, excluding unlabeled pixels.

The networks SegNet, PSPNet, DeepLabV3+, UNet and ResUnet are adopted as the most commonly used semantic segmentation networks in this study. SegNet is proposed by Badrinarayanan et al. [2017], and its network structure is a convolutional encoder-decoder. It mainly uses 2D convolution, 2D max pooling, and 2D max upsampling pooling. The input pooling indices of 2D max upsampling pooling come from the corresponding output of 2D max pooling. The input and output features of 2D convolution depend on the input and output data, respectively. The output feature number of the last convolution is the number of cloud-type categories. Under the function of softmax and argmax, the cloud type corresponding to the pixel can be obtained.

PSPNet is proposed by Zhao et al. [2017] and ResNet50 [He et al., 2016] is used as a feature map extractor for PSPNet in this study. The size of the feature map is reduced to 1/4 of the original input size to conserve memory. The feature maps in different levels are generated by the pyramid pooling module, and all feature maps are concatenated to predict cloud types.

DeepLabV3+ is proposed by Chen et al. [2018] and the feature map extractor of DeepLabV3+ is the same as that of PSPNet. In encoder processing, the encoder mainly uses four convolutional kernels with different dilation parameters and one global average pooling. In decoder processing, the low-level features obtained by the feature map extractor and the features obtained by the encoder are concatenated to predict cloud types.

Ronneberger et al. [2015] designed UNet. UNet is widely used in the field of imaging [Waldner & Diakogiannis, 2020, Yu et al., 2023, Yoo et al., 2022, Amini Amirkolaee et al., 2022, Pang et al., 2023], and it has significant similarities with SegNet. The main difference is that UNet concatenates the feature maps corresponding to the downsampling process during the upsampling process, and ultimately predicts cloud types. Zhang et al. [2018] developed ResUnet, which is built with residual units and has similar architecture to that of UNet.

In this study, a novel simple and effective deep learning-based network for cloud-type classification, named CldNet, is proposed and illustrated in Fig. 7. Depthwise separable convolution (DWConv) is applied into CldNet, which factorizes a standard convolution into a depthwise convolution [Chollet, 2017] followed by a pointwise convolution [Hua et al., 2018]. The difference between atrous DWConv and DWConv is the dilation of depthwise convolution. CldNet mainly consists of a DW-ASPP module and a DW-U module in Fig. 7. The DW-ASPP module and DW-U module are inspired by the atrous spatial pyramid pooling (ASPP) module of DeepLabV3+ and UNet, respectively. The DW-ASPP module uses four atrous depthwise separable convolutions with different dilations to capture useful multi-scale spatial context, and this can enhance the receptive field of the feature map. The DW-U module can extract feature maps at different levels. The fusion of the two extracted feature maps is beneficial for mining the inherent relationships between model input data and cloud types and achieving cloud-type recognition. More importantly, the two modules of CldNet have fewer parameters, which saves memory and makes them easier to deploy on edge devices.

The probability distribution of cloud-type prediction for each pixel is expected to be consistent with that of the cloud-type reference. Cross entropy can measure the difference information between two probability distributions. In the multi classification problem of deep learning, cross entropy can be replaced by negative logarithmic likelihood to calculate model loss. The loss update model parameters through backpropagation to achieve optimal performance of the model. Negative logarithmic likelihood loss (NLLLoss) is computed by Eq. 1 and Eq. 2, and the label of cloud type is calculated by Eq. 3.

where $N$ is the number of samples per batch; $C$ is the number of cloud-type labels; $Cld_{ref}^{prob}$ and $[label]$ is the probability and label of the reference cloud type, respectively; $Cld_{pre}^{prob}$ and $Cld_{pre}^{label}$ are the probability and label of the model’s predicted results, respectively; $x$ is the output result of the network structure; and the function $argmax$ returns the indices of the maximum value of all elements in the input vector.

This section will provide a detailed introduction about experimental setup. Satellite spectral data and cloud-type data from 2020 and 2021 are used as the training dataset, and those from 2022 are used as the test dataset. During the training process, the training dataset is randomly divided into five parts. Four of these parts are used for training, whose main purpose is to update the parameters of the model. The remaining part is used for validation, whose purpose is to prevent overfitting by early stopping of the model training [Dietterich, 1995] when the loss on the validation dataset no longer decreases after 10 consecutive epochs. The test dataset is used to evaluate the performance of the trained model.

All models in this study are trained on single NVIDIA GPU A100. In each training batch, single remote sensing image is segmented into $5\times 5$ slices and input into the model, and its dimension is $25\times 76\times 480\times 480$. The batch size is 25. Optimization algorithm Adam is used and its learning rate (lr) is controlled through a multi-step learning rate approach. If epoch $<$ 5, lr = 0.01; If 5 $\leq$ epoch $<$ 10, lr = 0.001; If 10 $\leq$ epoch $<$ 20, lr = 0.0001; If 20 $\leq$ epoch $<$ 30, lr = 0.00001; And if epoch $\geq$ 30, lr = 0.000001.

The recognition of nine cloud types and clear sky is a multi classification problem. In order to better evaluate the predictive classification performance of the model, accuracy is a commonly used indicator. $\mathrm{Accuracy}$ can measure the ratio of correctly classified predictions to the total number, as follows:

In order to understand the model’s ability to distinguish between clouds and clear skies, the indicator accuracy (N/Y for cloud) is defined as follows:

$\mathrm{Precision}$ focuses on evaluating the proportion of ture positive data among all predicted positive data. $\mathrm{Recall}$ focuses on evaluating how much data has been successfully predicted as positive among all positive data. For $\mathrm{precision}$ and $\mathrm{recall}$, each class needs to calculate its precision and recall separately. $\mathrm{F_{1}\mbox{-}score}$ is a comprehensive indicator of both $\mathrm{precision}$ and $\mathrm{recall}$, as follows:

$\mathrm{F_{1}\mbox{-}score_{macro}}$ directly adds up the $\mathrm{F_{1}\mbox{-}score}$ of different classes to calculate the average, which can treat each class equally.

$\mathrm{F_{1}\mbox{-}score_{weight}}$ is the sum of the $\mathrm{F_{1}\mbox{-}score}$ of each class multiplied by its weight, and this method considers class imbalance issues.

$\mathrm{F_{1}\mbox{-}score_{micro}}$ adds the TP, FP, and FN of each class first, and then calculates them based on the binary classification.

where TP is true positive; FN is false negative; FP is false positive; and TN is true negative. True and false refer to the correctness of the test result. True means that the test result is correct, and false means that the test result is incorrect. Positive and negative refer to the test results of the sample. Positive means that the intended target is detected, and negative means that the intended target is not detected. $N$ is the total number of samples in sigle image. $C$ is the number of cloud-type labels. $\mathrm{R}_{j}$ is the true distribution proportion of the class $j$.

In the previous section, our proposed model CldNet is found to perform best in all models, and thus CldNet is chosen for a more in-depth study. In order to improve the accuracy of the model, some auxiliary information has been added. The auxiliary information includes satellite zenith angle (SAZ), satellite azimuth angle (SAA), solar zenith angle (SOZ), solar azimuth angle (SOA), longitude, latitude, altitude, and underlying surface attributes (land or water). The impact of the auxiliary information on CldNet is explored in Table 2, and it is found that CldNet-W, which adds SAZ, SAA, SOZ and SOA to the model inputs, performs best. Compared with CldNet, the accuracy of CldNet-W has been improved by approximately 1.35%. This indicates that the auxiliary information can enhance the predictive ability of the model to identify cloud types.

An important aim of this study is to extrapolate cloud types in the nighttime areas by training model parameters using daytime label data. VIS and NIR channels can be used during the day, but not at night. If the model does not use VIS and NIR information, mask module will set the data involving VIS and NIR channels to zero. CldNet-O uses mask module to set the model input data involving VIS and NIR channels to zero through mask ratio and mask bands based on CldNet-W.

The overall cloud-type distributions predicted by the trained CldNet-W and CldNet-O at 2022-09-23 03:00 (UTC+0) are shown in Fig. 12, whose specific indicator results are recorded in Table A.2 and Table A.4. In Fig. 12(c), CldNet-W performs better than CldNet-O in cloud-type identification for most regions. However, CldNet-O extrapolates cloud types in the nighttime region due to not involving VIS and NIR channels. Compared with CldNet-W, the accuracy/$\mathrm{F_{1}\mbox{-}score_{micro}}$ and accuracy (N/Y for cloud) of CldNet-O have decreased by 7.34% and 0.45%, respectively. It can be seen that the removal of VIS and NIR has little effect on the distinction between clear and cloudy skies, but a significant impact on the identification of cloud types, and this result is also confirmed by the study of Zhang et al. [2019]. In order to construct a high-fidelity, all-day, spatiotemporal cloud-type database with spatial resolution $0.05^{\circ}\times 0.05^{\circ}$ over the entire satellite observation area, it is vitally important that CldNet-W and CldNet-O can be deployed online for daytime and nighttime areas, respectively.

The overall distributions of cloud types every 4 h for CldNet-O and reference during 2022-09-22 are shown in Fig. 13. Compared to the reference in Fig. 13(b), the cloud-type distributions in the nighttime area can be observed using CldNet-O in Fig. 13(a). In order to evidence the predicted results in the nighttime area, a method [Shang et al., 2017] for the identification of clear and cloudy sky over land based on BT of the satellites Himawari-8/9 spectral channel B14 is adopted. Here, confidence of clear sky for each pixel in the satellites Himawari-8/9 image is computed through the formula $\mathrm{(B14-270)/(288-270)\times 100\%}$. The overall distributions of confidence of clear sky every 4 h during 2022-09-22 are shown in Fig. 14.

The cloud cover over Australia (the black box in Fig. 13) is enlarged and compared with the confidence of clear sky at 2022-09-22 16:00 (UTC+0) in Fig. 15. The confidence of clear sky and cloud cover in the region R01 in Fig. 15 exhibit a very similar distribution, which verifies the model’s ability to capture cloud patterns. The cloud cover in the region R02 in Fig. 15(a) shows a long strip distribution, which is supported by the confidence of clear sky in Fig. 15(b). The cloud cover and confidence of clear sky in the region R03 are generally consistent. Overall, the results indicate that CldNet-O can capture subtle features of cloud distribution.

In this section, we directly apply the previously trained model, i.e., the trained CldNet-W, and CldNet-O without any fine-tuning to satellite spectral data with spatial resolution $\mathrm{0.02^{\circ}\times 0.02^{\circ}}$ from Himawari-8/9 satellite sensors to obtain the cloud-type distributions with the same spatial resolution $\mathrm{0.02^{\circ}\times 0.02^{\circ}}$. Four time points, including 2022-03-20 2:00 (UTC+0), 2022-06-21 2:00 (UTC+0), 2022-09-22 2:00 (UTC+0), and 2022-12-21 2:00 (UTC+0), are selected for cloud-type prediction. The cloud-type distributions with spatial resolution $\mathrm{0.02^{\circ}\times 0.02^{\circ}}$ predicted by the trained CldNet-W and CldNet-O are shown in Fig. 16(a) and 16(b), respectively. The JAXA’s P-Tree system only provides the cloud-type product with spatial resolution $\mathrm{0.05^{\circ}\times 0.05^{\circ}}$, which is used as a reference in Fig. 16(c). From the overall distribution of cloud types in Fig. 16, the trained CldNet-W and CldNet-O have achieved good performance.

In order to quantitatively evaluate the generalization ability of our models CldNet-W and CldNet-O, two methods, Low2High and High2Low, are used to calculate cloud-type classification indicators due to the resolution difference between prediction and reference. The process of Low2High is that the reference is first interpolated from $0.05^{\circ}\times 0.05^{\circ}$ to $0.02^{\circ}\times 0.02^{\circ}$ resolution, and then this is compared to the prediction to calculate those classification indicators. The process of High2Low is that the prediction is first interpolated by neighboring interpolation from $0.02^{\circ}\times 0.02^{\circ}$ to $0.05^{\circ}\times 0.05^{\circ}$ resolution, and then this is compared to the reference to calculate those classification indicators.

The results of the cloud-type classification indicators are recorded in Table 3. The indicators obtained by the method High2Low are better than those obtained by Low2High. For accuracy (N/Y for cloud), the results of CldNet-W and CldNet-O are all above 88%, which indicates that the model has excellent ability in the distinction between clear and cloudy skies. Overall, the trained models CldNet-W and CldNet-O have achieved good results in directly applying high-resolution satellite spectral information to predict high-resolution cloud-type distributions in term of accuracy/$\mathrm{F_{1}\mbox{-}score_{micro}}$, $\mathrm{F_{1}\mbox{-}score_{macro}}$, $\mathrm{F_{1}\mbox{-}score_{weight}}$, and accuracy (N/Y for cloud). This also demonstrates the good generalization ability of the model on higher resolution model input data.

One limitation of this study is that the satellite Himawari-8 was replaced by the satellite Himawari-9 in December 2022. Therefore, the training data is from the satellite Himawari-8, while the test data is partially from the satellite Himawari-9. Although the sensors carried by Himawari-8 and Himawari-9 satellites are the Advanced Himawari Imager, the spectral data distribution of the two satellites may be slightly different. The experimental results in Section 4.1.1 indicate that directly applying the model trained on satellite spectral data from the Himawari-8 sensor to satellite spectral data from the Himawari-9 sensor will reduce the accuracy of cloud-type recognition. In order to ensure accuracy, it is necessary to train the spectral data obtained from each satellite separately to obtain the model parameters of the corresponding satellite spectral data. Another limitation is that the labels used for model training are daytime cloud-type labels. In order to enhance the credibility of the model, model parameters still need to be trained using other data containing nighttime cloud-type labels.

In future research, we will adopt appropriate measures to optimize the limitations mentioned above and strive to expand the application scope of this model from regional to global scale. In order to achieve global cloud-type coverage, the proposed model will be applied to multiple geostationary observation satellites, such as Meteosat, GEOS-W/E, FengYun, and Himawari. The global distribution of cloud types is crucial in global climate change and environmental assessment research.