A Machine Learning-Based Quality Control Algorithm for Heavy Rainfall Using Multi-Source Data

1. Introduction

In surface meteorological observation, precipitation data are regarded as one of the most frequently utilized and critically supportive datasets [1,2,3]. Accurate precipitation data, particularly heavy rainfall data, is widely employed in weather forecasting, climate change studies, and numerical model assimilation [3,4,5]. These data also serve as essential decision-making bases for key societal sectors such as agricultural production, water resource management, and ecological security [6,7,8].

With the continuous construction and expansion of Chinese surface observation network, the number of automatic weather stations (AWSs) has exceeded 70,000, significantly improving observational coverage and density [9]. However, more than 95% of them are unmanned. Without dedicated on-site maintenance or backup precipitation observation devices, these unmanned stations often suffer from operational instability and non-standard installation, making them prone to various data quality anomalies [10,11,12].

Currently, the quality control (QC) of surface precipitation data in China primarily relies on traditional methods, such as regional extreme value checks, temporal consistency checks, and spatial consistency comparisons between stations [13,14,15,16]. For instance, the Meteorological Data Operation System (MDOS) of the China Meteorological Administration (CMA) employs a combination of file-level rapid QC, hourly batch QC, and daily data QC to achieve automated processing [17]. However, these methods still exhibit high rates of missed detections and false alarms in cases such as isolated heavy rainfall caused by human operational errors or continuous large precipitation values resulting from equipment malfunctions [18,19,20,21,22]. Although manual review can serve as an effective supplement to automated QC, it is typically conducted several hours after data arrival, making it difficult to meet the timeliness requirements for real-time heavy rainfall monitoring and emergency response [17,23].

With the advancement of ground- and satellite-based observation technologies and the explosive growth of meteorological data, multi-source data collaborative QC for precipitation has demonstrated significant advantages [24,25,26,27,28]. Compared with traditional approaches, these methods integrate various observational data, including Doppler weather radar, geostationary meteorological satellite, and weather detector, effectively reducing the false alarm rate in anomaly detection [24,25,26]. Furthermore, multi-source collaborative QC also contributes positively to the quality improvement of hydrological precipitation data [27,28]. However, multi-source collaborative QC methods rely mainly on statistically based static thresholds and have not yet fully exploited the complementary potential of multi-source data. Therefore, their optimization and fusion mechanisms require further exploration.

Machine learning methods can automatically learn discriminative features from historical data and identify complex patterns that are difficult to capture using static thresholds [21,29,30]. This capability enables the establishment of more robust and complex nonlinear mapping relationships between precipitation data and other multi-source observations. Additionally, through model structure design or ensemble learning strategies, machine learning can adaptively adjust the contribution weights of various inputs, effectively integrating the strengths of different data sources [31,32]. For example, machine learning algorithms such as decision trees, K-nearest neighbors, and isolation forests have been applied to classify outliers, significantly improving the detection accuracy of anomalous precipitation data [33,34,35,36]. Neural networks have been used to construct real-time precipitation confidence intervals, enabling model adaptation across spatial and seasonal variations [37,38]. By incorporating topographic elevation data through deep neural networks and multi-scale ensemble learning techniques, the capability to detect precipitation quality anomalies in data-sparse regions such as coastal and mountainous areas has been enhanced [39].

To address the missed detections and false alarms in the heavy rainfall QC caused by human operational errors, equipment malfunctions, and other non-meteorological factors, a machine learning-based algorithm was developed in this study. The proposed algorithm employed multi-source data (AWS, radar, satellite, and metadata) from 1 June 2022 to 31 December 2024 and was modeled using advanced machine learning techniques. Furthermore, interpretability analysis was performed on the model results to enhance their reliability, thereby providing more dependable decision support for weather forecasting operations and related public welfare sectors.

2. Data and Methods

2.1. Data

The data employed in this study are categorized into three primary types: (1) surface observational data (including hour-level and minute-level data), (2) remote sensing data, and (3) metadata. The dataset covers a two-year period from 1 June 2022 to 31 May 2024 for model training and testing, and a subsequent six-month period from 1 June to 31 December 2024 for validation, with the spatial scope encompassing the entirety of China. For the development of a QC algorithm for heavy rainfall, this study focused on surface observational data with hourly precipitation exceeding 40 mm. The corresponding remote sensing data and metadata from these events were selectively extracted for analysis, with a total dataset of over 127,000 records.

2.1.1. Surface Observational Data

The surface observational data were obtained from a CMA-maintained network comprising national and provincial AWSs [6,9,17], with a station density higher in the eastern and southern regions than in the western and northern areas (Figure 1). Thirteen hour-level parameters, including hourly precipitation (PRE_1h), temperature (TEM), dew point temperature (DPT), atmospheric pressure (PRS), 3-hour atmospheric pressure change (PRS_Change_3h), 24-hour atmospheric pressure change (PRS_Change_24h), relative humidity (RHU), instantaneous wind speed (WIN_S_INST), maximum instantaneous wind speed in 1 hour (WIN_S_INST_Max), 10-minute average wind speed (WIN_S_Avg_10mi), instantaneous wind direction (WIN_D_INST), direction of maximum instantaneous wind speed in 1 hour (WIN_D_INST_Max), and 10-minute average wind direction (WIN_D_Avg_10mi), were extracted as input features for the models.

In addition to the hour-level surface observational data, minute-level precipitation data were also incorporated in this study. Due to the substantial volume of the minute-level data, feature extraction processes were performed to characterize its mean, extremums, and degree of dispersion. The average (PRE_1m_Ave), maximum (PRE_1m_Max), minimum (PRE_1m_Min), and standard deviation (PRE_1m_Std) of all minute-level values in each hour were calculated and used as input features for the models. The average and standard deviation of minute-level precipitation were defined by Equations (1 and 2), respectively.

$$PRE\_1m\_Ave={\displaystyle \frac{1}{N}}{\sum }_{i=1}^N{PRE}_i$$

(1)

$$PRE\_1m\_Std=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{i=1}^N{({PRE}_i-PRE\_1m\_Ave)}^2}$$

(2)

where ${PRE}_i$ represents the minute-level precipitation value of the i-th sample, and N indicates the record number of minute-level precipitation in this hour.

2.1.2. Remote Sensing Data

The remote sensing data utilized in this study comprised the composite reflectivity (CR) mapping product from Doppler weather radars and the cloud-top temperature (CTT) product retrieved from the Fengyun-4B (FY-4B) geostationary meteorological satellite. Both products were quality-controlled data [40,41].

The CR mapping product is generated through periodic volume scans conducted by individual radar stations [42]. Base data, such as reflectivity, are transmitted in real time to processing centers, where they undergo rigorous QC algorithms to remove ground clutter and electromagnetic interference. The polar coordinate data from each single radar are then interpolated onto a unified latitude–longitude grid. Finally, a “maximum reflectivity“ fusion algorithm is applied to produce the composite reflectivity mapping product. This mapping product is built upon a foundation of unified calibration across all radars in China, ensuring consistency and comparability of data from different stations. It provides complete coverage of China at high spatiotemporal resolutions of 6 minutes and 0.01°.

The FY-4B satellite-derived cloud-top temperature product is retrieved via multi-channel remote sensing detection [43]. The satellite employs an Advanced Geosynchronous Radiation Imager to observe specific infrared channels, such as the long-wave infrared window channels at 10.3-11.3 μm and 11.5-12.5 μm. The raw radiance emitted from cloud tops is acquired and radiometrically calibrated to obtain accurate top-of-atmosphere radiance values [44]. Based on Planck’s blackbody radiation law, these radiance values are subsequently inverted into cloud-top temperatures. The geographic projection coordinates from the satellite data are computed based on the World Geodetic System-1984 Coordinate System [45], and the row and column numbers of the detection products are thereby converted into geographical latitude and longitude values. This product covers East Asia with temporal and spatial resolutions of 15 minutes and ~4 km, respectively.

Radar and satellite data exhibit strong complementary characteristics. Radar observations provide high spatiotemporal resolution, enabling effective monitoring of the initiation and evolution of meso- and micro-scale severe convective systems, thereby compensating for the limitations of satellites in detailed monitoring. Satellite observations are not constrained by topography and can cover regions such as oceans, plateaus, and mountains, where radar deployment is challenging, thus mitigating radar blind zones in areas with complex topography. By integrating radar and satellite data, the technical advantages of both ground-based and satellite-based remote sensing can be leveraged, enhancing the reliability of heavy rainfall QC.

2.1.3. Metadata

Metadata are defined as “data about data”, which describe the context and provenance of the data and are essential for their correct interpretation and use. Key components include station location, instrument specifications, measurement method, etc. [46]

The acquisition of metadata is a systematic process. During station construction, geographical and instrumental information are surveyed and recorded. Subsequently, identifiers and timestamps are embedded by the data loggers during data collection and transmission. Finally, processing logs are integrated through standardized formats in processing and archiving. In this study, altitude (Alti), station level, and measurement method were selected as the input features for the models.

2.2. Feature Engineering

2.2.1. Data Matching

Hourly precipitation data represent the cumulative rainfall during the one-hour period immediately preceding each hour (e.g., the precipitation recorded at 06:00 corresponds to the accumulated rainfall between 05:00 and 06:00). However, remote sensing data sources—with temporal resolutions of 6 minutes for radar and 15 minutes for satellite—exhibit significantly higher temporal resolution than hourly precipitation data. This discrepancy in temporal scale leads to fundamental inconsistencies when directly performing grid-to-grid matching in these heterogeneous datasets.

A sliding average temporal window method was applied to process the remote sensing data [47]. For each hourly precipitation value, all available radar or satellite gridded observations within the preceding one-hour window were extracted. The average of these high-resolution values was then computed (A schematic diagram is shown in Figure 2) and matched with the hourly precipitation data. This method effectively aligns the high-frequency remote sensing data with the hourly precipitation observations on a consistent temporal scale, thus achieving temporal consistency for multi-source data matching.

For spatial matching, a spatial nearest-neighbor averaging method was employed in this study [48]. Centered on the latitude and longitude coordinates of each AWS, the CR and CTT gridded data within the surrounding area were searched. The arithmetic mean of the five nearest grid points was then calculated and used as the remote sensing observation values matched to the hourly precipitation data.

The utilization of multiple grid points within a defined spatial domain for averaging is motivated by practical considerations regarding spatial representation and error minimization. During heavy rainfall events, accompanying strong low-level winds can induce horizontal advection of hydrometeors, resulting in a displacement between the precipitation measured at ground stations and the cloud properties observed directly overhead via remote sensing. Relying solely on the single grid value immediately above one AWS may introduce significant mismatches. This method reduced small-scale spatial inconsistencies, thereby providing a more reliable representation of the remote sensing data within the vicinity of the station.

2.2.2. Feature Selection

In high-dimensional datasets, redundancy and multicollinearity among input features are frequently observed [49]. These issues can increase computational complexity, prolong training time, and potentially lead to overfitting, thus reducing model generalization capability and interpretability. Therefore, feature selection is a critical preprocessing step in machine learning. It aims to construct a low-dimensional and efficient feature subset while preserving essential information from the original data. This process enhances model performance and computational efficiency.

This study adopted a strategy based on correlation analysis to systematically eliminate redundant features [50]. The Pearson correlation coefficient matrix was calculated for all pairwise combinations of the input features and a threshold was set to identify highly correlated feature pairs. For each such pair, one of the features was retained while the other was eliminated. This process effectively removed redundant information, ensuring that the final feature subset consists of features with low mutual correlation.

Based on absolute Pearson correlation coefficient analysis (Figure 3) with a threshold of |r| > 0.5, a systematic elimination of redundant features from the initial 22 input features were conducted. This process resulted in the removal of 7 redundant features (Table 1), including altitude (Alti), dew point temperature (DPT), 24-hour pressure change (PRS_Change_24h), instantaneous wind speed (WIN_S_INST), instantaneous wind direction (WIN_D_INST), minute-level precipitation mean (PRE_1m_Ave) and standard deviation (PRE_1m_Std) in an hour . The refined feature set retained 15 core features, achieving a 31.8% reduction in feature dimensionality while fully preserving critical precipitation-related information encompassing dynamic, thermodynamic, moisture, and cloud microphysical characteristics.

2.3. Machine Learning Models

Gradient boosting machine (GBM) is an ensemble learning framework that leverages the collective power of multiple base learners (typically decision trees) [51,52]. It works by iteratively training a sequence of these base learners to enhance prediction performance. This framework is capable of effectively capturing complex non-linear relationships among variables, adapting to datasets of various structures and sizes, and is widely applicable to both regression and classification tasks.

Based on GBM framework, four particularly prominent and widely-used models have emerged: eXtreme Gradient Boosting (XGBoost) [53,54], Light Gradient Boosting Machine (LightGBM) [55], Categorical Boosting (CatBoost) [56,57], and Gradient Boosted Regression Trees (GBRT) [58]. These models exhibit distinct advantages in areas such as training efficiency, handling of missing values, and mitigation of overfitting. In this study, these four advanced models were employed to systematically train and predict on the dataset. To ensure these models achieve their designed performance, hyperparameter optimization was conducted automatically using the Optuna library [59] in Python. This library supports efficient hyperparameter tuning across multiple machine learning models and was utilized to identify optimal parameter configurations.

Since the input data are derived from multiple heterogeneous sources, significant differences are observed in their measurement principles, accuracy, and acquisition environments. Considerable variations are exhibited in the value ranges, distribution characteristics, and systematic biases among different inputs [60]. Therefore, Z-score normalization was applied to preprocess each type of input feature separately, by using the formula given in Equation (3).

$$z_i^j={\displaystyle \frac{X_i^j-{\mu }_i}{{\sigma }_i}}$$

(3)

where $z_i^j$ denotes the Z-score of the j-th sample in the i-th feature, $X_i^j$ represents the observed value of the the j-th sample in the i-th feature, ${\mu }_i$ is the mean value of the i-th feature, and ${\sigma }_i$ is the standard deviation of the i-th feature.

By subtracting the mean and dividing by the standard deviation, the features were transformed into a unified scale, thereby effectively reducing distribution discrepancies. Z-score normalization prevented certain features from exerting excessive influence on model training due to their large numerical values, while also contributing to accelerated convergence and improved generalization capability.

2.4. Model Performance Metrics

To evaluate model performance, a confusion matrix-based performance metrics employed in this study (Figure 4), utilizing accuracy, precision, recall, and F1-Score as key metrics. These metrics are defined by Equations (4–7), respectively [61,62].

$$\mathrm{a}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{y}={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(4)

$$\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}={\displaystyle \frac{TP}{TP+FP}}$$

(5)

$$\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{l}={\displaystyle \frac{TP}{TP+FN}}$$

(6)

$$F_1\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=2\times {\displaystyle \frac{Pre\times Rec}{Pre+Rec}}$$

(7)

where TP (true positive) denotes the number of samples that are positive and correctly predicted as positive; FN (false negative) refers to the number of samples that are positive but incorrectly predicted as negative; FP (false positive) is the number of samples that are negative but incorrectly predicted as positive; and TN (true negative) represents the number of samples that are negative and correctly predicted as negative.

The “actual value“ in the confusion matrix is obtained from manually corrected hourly precipitation QC results, which are currently regarded as the most accurate benchmark (although its timeliness is relatively poor). Furthermore, by comparing the metrics calculated between the model’s outputs and the actual value against those derived from traditional MDOS results and the same actual value, the proposed GBM machine learning models effectively demonstrate their superiority over current operational QC methods.

2.5. Model Interpretation

SHapley Additive exPlanations (SHAP) is widely regarded as one of the most highly useful tools in model interpretability [63,64]. It is a method founded on cooperative game theory that is used to interpret the predictions generated by machine learning models. The prediction of the black-box model is decomposed into the sum of individual feature effects by calculating each feature’s contribution, and renders the final prediction clear and interpretable. A lower SHAP value indicates a lower feature contribution, whereas a higher value corresponds to a greater contribution. The SHAP value is mathematically defined by Equation (8).

$$g(x^{\prime })={\varnothing }_0+{\sum }_{i=1}^M{\varnothing }_i{x}_i^{\prime }$$

(8)

where $g\left(x^{\prime }\right)$ represents the explanation model, x denotes the input features, and M indicates the number of input features. ${\varnothing }_i$ represents the Shapely values for feature i and ${\varnothing }_0$ denotes the constant output value when all inputs are absent. For each feature, the term ${\varnothing }_i$ is computed using Equations (9) and (10).

$$(f,x)={\sum }_{z^{\prime }\subseteq x^{\prime }}{\displaystyle \frac{\left|z^{\prime }\right|!(M-|z^{\prime }|-1)!}{M!}}[f_x(z^{\prime })-f_x(z^{\prime }\backslash i)]$$

(9)

$$f_x\left(z^{\prime }\right)=E\left[f\right(x\left)\right|x_S]$$

(10)

where f represents the black box model. $z^{\prime }$ denotes a subset of input features, and $\left|z^{\prime }\right|$ is the number of non-zero entires in $z^{\prime }$. $E\left[f\right(x\left)\right|x_S]$ represents the expected value of the function and S is the set of non-zero indices in $z^{\prime }$.

SHAP values can be approximated using various methods, such as Kernel SHAP, Deep SHAP, and Tree SHAP. In this study, the Tree SHAP-based approach was employed and implemented via the SHAP library in Python.

3. Results

3.1. Importance Analysis

The analysis of normalized feature importance in the four models for the 15 features revealed that remote sensing data played a dominant role in the predictions, with mean normalized feature importance values of 0.332 for CR and 0.220 for CTT, respectively (Figure 5). Precipitation-related features were consistently identified as stable predictors, where the PRE_1m_Max, PRE_1h, and PRE_1m_Min were identified as subdominant predictive factors, with mean normalized feature importance values of 0.090, 0.064, and 0.056. Moderate contributions were observed from wind speed, wind direction, and TEM, with mean normalized feature importance values ranging from 0.026 to 0.046. In contrast, the importance of RHU, PRS, PRS_Change_3h, and metadata features were relatively low, with normalized feature importance values all below 0.01.

However, distinct model-specific characteristics were observed in the feature weighting in the four models (Figure 6). In the XGBoost model, the ranking of feature importance was observed to be similar to the average of the four models, with remote sensing data identified as the most significant (CR, 0.363; CTT, 0.232), followed by precipitation-related features (PRE_1m_Max, 0.080; PRE_1m_Min, 0.067; PRE_1h, 0.056). Greater sensitivity to wind direction (WIN_S_Avg_10mi, 0.113; WIN_D_INST_Max, 0.051) and wind speed (WIN_S_INST_Max, 0.063; WIN_S_Avg_10mi, 0.056) features were exhibited by the LightGBM model, with these values being significantly elevated compared to other models. A unique pattern was presented by CatBoost, which was identified as the only model where a higher importance was allocated to CTT (0.313) than to CR (0.257). This model was also characterized by a considerably greater emphasis being placed on the WIN_S_INST_Max (0.078). In contrast, an extreme dependence on CR was displayed by the GBRT model, which was attributed a normalized feature importance value of 0.509, accounting for over half of the total feature importance. Concurrently, wind-related features were largely disregarded (normalized feature importance value≤0.006). These discrepancies highlighted the inherent biases of different models in capturing complex feature-target relationships, underscoring that the selection of an appropriate model is critical for the accurate interpretation of feature influences in practical applications.

During the modeling process, feature elimination was performed based on feature importance. Only features whose cumulative contribution exceeded 99.5% after sorting in descending order were retained. This procedure further streamlined the model architecture and improved interpretability, following the initial feature selection outlined in Section 2.2.2.

3.2. Comparison of Performance in Different Models

The performances of the four models in heavy rainfall QC were systematically evaluated by the metrics outlined in Section 2.4. The original dataset was partitioned into training and testing sets at a 7:3 ratio. Comparative analysis of two sets revealed that the performance degradation of the four models was not significant, with observed decreases in F1-scores of 0.045, 0.049, 0.047, and 0.044 respectively, thus indicating no substantial overfitting were detected (Table 2 and Table 3).

High accuracies (accuracy > 0.99) were observed by all models in testing set (Table 3). However, it could be attributed to the significant class imbalance in the precipitation dataset employed, where normal precipitation events (negative samples) substantially outnumbered anomalous precipitation events (positive samples). The identification of true negative samples by the models was relatively easy, resulting in the TN value in Eq. (2) being considerably larger than TP, FP, and FN. Consequently, the accuracies were universally inflated, failing to adequately reflect the capability to detect the critical minority class of anomalous precipitation. In comparison, the precision, recall, and F1-score demonstrated greater diagnostic value in identifying such imbalanced datasets.

According to the results, all evaluated metrics of the four models outperformed those of the MDOS QC (Table 3). The superiority of applying machine learning techniques to multi-source data, particularly remote sensing data, for heavy rainfall QC was effectively highlighted. The highest F1-score (0.816) and recall (0.789) were achieved by XGBoost, indicating its superior comprehensive performance in identifying anomalous precipitation. LightGBM was observed to closely follow, with its F1-score (0.811) nearly matching that of XGBoost, while a slightly higher precision (0.845) was attained. The highest precision (0.849) was recorded by the GBRT model, though its recall was found to be relatively lower. All metrics obtained by the CatBoost model on the test set were noted to be slightly inferior to those of the other three comparative models.

The precision-recall (PR) and receiver operating characteristic (ROC) curves [65] of the four models were observed to be similar (Figure 7). The highest average precisions (APs) of 0.777 were achieved by both XGBoost and CatBoost. Meanwhile, the highest two area under the curves (AUCs) were exhibited by LightGBM and XGBoost, with values of 0.942 and 0.941, respectively. Therefore, through a synthesis of metrics, RP and ROC curves assessment, the XGBoost model was regarded as the preferred high-quality model and was chosen for further study (Section 3.4). Compared to the MDOS, the XGBoost model achieved an increase in precision by 0.110, recall by 0.162, and F1-score by 0.140 in the heavy rainfall QC.

3.3. Model Validation

To evaluate the actual predictive performance of the models, data from 1 June to 31 December 2024 were used as the validation set, with the performance metrics calculated accordingly (Table 4). It was observed that slight decreases in precisions were exhibited by all models on the validation set. Reductions of 0.006, 0.022, 0.010, and 0.025 were recorded for XGBoost, LightGBM, CatBoost, and GBRT, respectively, with more pronounced declines being shown by LightGBM and GBRT. A more substantial overall declines were demonstrated in recalls, with decreases of 0.062, 0.053, 0.066, and 0.056, respectively, and the largest reduction was observed in CatBoost. Influenced by both precisions and recalls, a consistent downward trends in the F1-Score were also presented on the validation set, with declines ranging between 0.037 and 0.043 being exhibited by the four models. The overall degradation in model performance may be attributed to the temporal distribution discrepancy between the validation and training sets, indicating a limited decline in generalization capability when models are confronted with temporally shifted data. Nevertheless, all models are found to significantly outperform the MDOS QC results across all metrics in validation set, suggesting that satisfactory QC capability and practical utility are still maintained in actual prediction.

3.4. Feature Contribution

Rooted in an individualized model interpretation method, SHAP values allow a distinct explanation to be provided for each feature. Figure 8 presents how a feature modulates its own contribution to the model output. In this figure, the color of the dot represents the value of that feature for the review (red: high, blue: low), and the position of the dot is the contribution of the feature on the review helpfulness. Positive SHAP values signify that the feature’s value contributes to the model’s prediction of “anomalous precipitation“, and a larger value indicates a stronger contribution. Conversely, negative SHAP values signify a contribution to the prediction of “normal precipitation“, with a smaller (more negative) value corresponding to a stronger contribution.

CR was regarded as the most discriminative feature in the model (Figure 6a and Figure 8). The SHAP summary plot demonstrated that lower CR values were significantly positively correlated with the probability of a record being classified by the model as “abnormal precipitation“ (Figure 8). Additionally, CR SHAP values exhibited a dependence on CTT and PRE_1h (Figure 9a,d). The samples with CR values below 20 dBZ, combined with either hourly precipitation greater than 70 mm or CTT exceeding 260 K, were strongly classified as “anomalous precipitation” by the model. These phenomena were highly consistent with the physical mechanisms governing heavy rainfall formation. Weak radar echoes imply underdeveloped convective systems and weaker processes of water vapor condensation and hydrometeor growth, which are theoretically insufficient to support precipitation rate greater than 40 mm/h. Consequently, when an observed hourly precipitation with a lower CR, the record is highly likely to be attributed to manual operational errors, such as garden watering, cleaning the gauge without disconnecting communication [20,25,26].

CTT was established as a key indicator in the model (Figure 6a and Figure 8). It showed that samples with elevated CTT (i.e., warmer cloud tops), particularly when combined with either a CR below 20 dBZ or hourly precipitation exceeding 70 mm, were readily classified as anomalous precipitation (Figure 10a,d). From a meteorological perspective, cloud systems capable of generating intense convection are typically characterized by high cloud top heights and low cloud top temperatures. A heavy rainfall event accompanied by a high CTT value is considered suspicious. Such events may be associated with shallow or warm-cloud-dominated systems, which typically lack the dynamic and microphysical conditions necessary for producing substantial rainfall. Similar to the CR SHAP analysis, these phenomena may also be attributed to non-meteorological factors.

The PRE_1m_Max was identified as a critical feature utilized for the model (Figure 6a and Figure 8). SHAP analysis indicated that the model was more likely to classify a precipitation event as anomalous under two scenarios: (1) the PRE_1m_Max was less than 0.5 mm (Figure 8 and Figure 11). This was primarily attributed to a logical inconsistency in the data, as such a low precipitation rate was physically implausible to support the observed hourly precipitation greater than 40 mm. This discrepancy is likely caused by data encoding or transmission errors [22,26]; (2) the PRE_1m_Max exceeded 3.5 mm, particularly when accompanied by either a CR of less than 20 dBZ or a CTT greater than 260 K (Figure 11a,b). An actual heavy rainfall event is composed of several high-intensity minute-scale segments, yet its maximum value usually falls within a reasonable physical range. This anomalously high minute values were likely to be induced by instrumental malfunctions [11,18,22,26], such as the persistent false triggering of a reed switch in a tipping-bucket rain gauge due to metal fatigue or foreign object adhesion, leading to an abnormally high count within a single minute. Additionally, occasional strong electromagnetic interference can also produce such physically implausible extreme peaks [20].

The PRE_1m_Min played a distinctive role in this model (Figure 6a and Figure 8). The SHAP analysis revealed that PRE_1m_Min did not heavily depend on other features (Figure 12). But anomalously low values (particularly negative values) of this feature served as a strong indicator for the model to classify a record as “abnormal precipitation”. As minute precipitation represents an accumulated quantity over a period, its value is inherently non-negative. Therefore, the occurrence of negative values exposes errors in the data encoding or transmission, such as improper sign bit handling or data type conversion errors [22,26]. Thus, it is also established as a sensitive indicator for heavy rainfall QC.

4. Summary and Conclusion

In this study, a machine learning-based QC algorithm for heavy rainfall was developed by integrating AWS observations with radar, satellite, and metadata inputs. Based on heavy rainfall samples from 1 June 2022 to 31 December 2024, performance metrics and interpretability were analyzed on the model outputs. The main conclusions are summarized as follows:

In heavy rainfall QC, the four gradient boosting models (XGBoost, LightGBM, CatBoost, and GBRT) demonstrated a marked improvement over the conventional MDOS method, achieving consistently superior performance in precision, recall, and F1-score. Through a comprehensive evaluation of the performance metrics, PR and ROC curves, the XGBoost model achieved the best overall performance, with precision, recall, and F1-score reaching 0.844, 0.789, and 0.816 on the testing set. Compared to the MDOS, the XGBoost model achieved an increase in precision by 0.110, recall by 0.162, and F1-score by 0.140 in heavy rainfall QC. This finding indicated that gradient boosting-based machine learning algorithms can effectively learn the nonlinear relationships between normal and abnormal precipitation from complex multi-source data. These algorithms significantly improved the identification of challenging QC issues like isolated anomalies and continuous heavy rainfall, reducing both false alarms and missed detections.

The radar, satellite, and minute-level precipitation data were found to play dominant roles in the model’s feature importance. The inclusion of these high-value features overcame the limitations of traditional QC methods that rely on static thresholds and single data sources. Moreover, radar data compensated for the spatiotemporal resolution limitations of satellite observations, while satellite data effectively covered the observational gaps of radar in complex terrain. Through the adaptive learning strategy of the machine learning model, these two data sources achieved complementary advantages and synergistic effects.

Interpretability analysis based on the SHAP method demonstrated that the model’s decision logic was highly consistent with meteorological physical mechanisms. For example, low CR or high CTT generally indicated weak precipitation or underdeveloped convective clouds. Heavy rainfall records associated with such features were accurately identified as anomalous by the model, which corresponds to manual operational errors (e.g., garden watering, cleaning the gauge without disconnecting communication) or other non-meteorological factors. Furthermore, anomalous minute-level precipitation extremes were confirmed by the model as key indicators for identifying instrument malfunctions (e.g. the persistent false triggering of a reed switch in a tipping-bucket rain gauge due to metal fatigue or foreign object adhesion), data encoding and transmission errors. This alignment between the decision logic mechanism and physical principles enhances the reliability of the model results and makes the decision process transparent and interpretable for forecasters, demonstrating strong potential for operational application.