This study focuses on Anhui Province in eastern China (116.5°–119.5° E, 29.5°–31.2° N) (Figure 1a), covering approximately 140,000 km². Located within the core region of the East Asian monsoon, the province experiences an average annual precipitation exceeding 1,200 mm, with pronounced spatiotemporal variability—over 70% of the annual total occurs from May to September (Figure 1c) (Wang et al., 2024; Wu et al., 2022). As a transitional zone between northern and southern climate regimes, Anhui is influenced by both cold air masses from Siberia and warm, moist flows from the Pacific Ocean, resulting in a complex atmospheric circulation system. This climatic configuration fosters frequent severe convective weather, including STHP events with intensities reaching 80–100 mm per hour (Hao et al., 2012).

The topography of Anhui is markedly heterogeneous (Figure 1b). The Yangtze and Huaihe Rivers divide the province into three distinct subregions: the Huaibei Plain in the north, characterized by flat terrain and agricultural landscapes; the central Jianghuai Hilly Area, which includes the western Dabie Mountains; and the southern Jiangnan Region, dominated by the mountainous terrain of Southern Anhui. Orographic lifting in the Southern Anhui and Dabie Mountains significantly enhances local precipitation intensity and frequency (Kim et al., 2022; Tong et al., 2017), often inducing secondary hazards such as flash floods and landslides. These precipitation patterns also contribute considerable uncertainty to flood risks in the Yangtze and Huaihe River basins.

The interplay of complex geography and monsoonal climate supports diverse dynamic–thermal coupling mechanisms for precipitation formation (Figure 1c). Combined with the high density of meteorological stations across the region (Figure 1b), Anhui Province offers an ideal natural laboratory for studying multi-scale physical processes of STHP in a climatic transition zone, thereby supporting improved early warning systems and regional disaster prevention strategies.

DNN is designed to extract hierarchical feature representations from raw inputs through successive nonlinear transformations. A standard layer-wise mapping can be summarized as:

In Equation 1, where x represents the input features, W and b denote the trainable weight matrix and bias vector, respectively, and σ is the nonlinear activation function. While powerful, standard DNN often encounter performance plateaus when processing multi-source meteorological data. In such scenarios, directly concatenating variables with disparate physical meanings and statistical scales (e.g., moisture vs. dynamic factors) can dilute the expressive capacity of individual feature streams and hinder the model’s ability to capture category-specific nuances.

To address these limitations, this study introduces a MSDNN framework (Figure 3). The core architectural philosophy is to decouple the initial feature extraction into five specialized sub-networks: SSI (severe convection and storm indices), HTV (height and thickness variables), HRV (humidity-related variables), SEI (stability and energy indices), and WVV (wind and vertical motion variables). Formally, the transformation for each sub-network i is expressed as:

In Equation 2, where f i denotes the nonlinear mapping function of the i -th sub-network, θ i represents its trainable parameters, and h i corresponds to the latent feature representation of that specific stream. These multi-stream features are subsequently integrated within a fusion layer, typically through concatenation or weighted combination:

In Equation 3, where g denotes the fusion function and θ fusion represents the parameters of the fusion stage. The resulting integrated representation is then fed into a fully connected layer, with the final identification result y generated through activation function mapping:

In Equation 4, where σ denotes the nonlinear activation function W represents the learnable weight matrix that linearly projects the fused feature vector and b is the trainable bias vector that shifts the activation threshold.

As illustrated in Figure 3, each sub-network channel is composed of stacked dense layers interspersed with batch normalization (BN) and Swish activation functions [ Swish ( x ) = x ∗ sigmoid ( x ) ] . The notation Dense ( d in → d out ) signifies a fully connected layer projecting an input dimension of d in to an output dimension of d out . Notably, the WVV stream is assigned a deeper initial architecture (two layers: 25 → 64 → 128) compared to other streams. This specific design is motivated by the high dimensionality (25 features) and the complex, multi-scale turbulent nature of wind field and vorticity data, which requires a larger hypothesis space to capture the non-linear fluid dynamics essential for STHP identification. Following three stages of stream-specific abstraction, the latent representations are integrated into a 1,152-dimensional fused vector, which is processed through two deep transformation layers (1,152 → 512 → 512) to model complex physical interactions. To ensure numerical stability and generalization, we implemented a multi-tier regularization strategy: (1) BN (Ioffe and Szegedy, 2015) is applied after each major transformation stage to stabilize training. (2) A Dropout layer (rate = 0.5) is inserted before the final output (Baldi and Sadowski, 2013), randomly deactivating 50% of neurons to prevent co-adaptation and enhance generalization. (3) The final layer projects the representation to a single unit (512 → 1) with a Sigmoid activation function to generate the probability of STHP occurrence.

The MSDNN was trained using the binary cross-entropy loss function and the Adam optimizer (Kingma and Ba, 2014). To facilitate convergence, an adaptive learning rate strategy (ranging from 10⁻⁵ to 10⁻³) was employed over a maximum of 10,000 epochs with a batch size of 32. To evaluate its performance, two ensemble learning algorithms—RF and LGB—were employed as benchmarks. RF utilizes a Bagging framework to reduce variance through collective decision-tree voting, while LGB implements a gradient boosting framework optimized by gradient-based one-side sampling (GOSS) and exclusive feature bundling (EFB) to enhance efficiency and handle large-scale datasets. These models serve as relevant baselines for assessing the MSDNN’s capacity to capture the complex nonlinear relationships inherent in STHP identification. Detailed hyperparameter specifications are provided in Table 2.

To quantify the marginal contribution of individual environmental physical variables to STHP discrimination, we employed the IG method for feature attribution analysis in the MSDNN model. Introduced by Sundararajan et al. (2017), IG satisfies fundamental axiomatic properties including sensitivity and implementation invariance, while offering completeness in attribution. This approach overcomes limitations of gradient-based methods in saturated regions through path integration.

Assume that the output probability of a MSDNN for the positive “STHP” class is F ( x ) ∈ [ 0 , 1 ] , where x i , i = 1 , 2 , 3 … n denotes the sample input and x i ′ represents the baseline value. IG is defined as the integral of the gradient of the output with respect to each feature along the straight-line path from the baseline to the sample:

In Equation 5, where IG i represents the attribution of the i -th feature to the sample x i , ∂ F ∂ x represents the gradient of the model output with respect to the input features. Using Riemann summation for approximation, m equidistant steps are sampled within the interval [ 0 , 1 ] , yielding Equation 6:

Owing to its completeness property, the sum of attributions across all features equals the total contribution to the model output F ( x ) relative to the baseline. This allows the model’s output probability to be decomposed into the baseline probability and the sum of individual feature attributions, yielding Equation 7:

In this study, the five-branch MSDNN processes 71 environmental variables for binary classification. IG computations were performed in the original input space following tensor preprocessing. For each sample, gradients were accumulated along integration paths from baselines to actual inputs across all five branches, generating feature-level attributions subsequently aggregated into sample-level spectra. To ensure numerical stability, stochastic components (e.g., Batch Normalization and Dropout) were disabled during inference. The integration step parameter m was set to 200, optimizing between computational efficiency and approximation accuracy.

The selection of baseline inputs critically influences IG attribution, as it establishes a reference state representing neutral or non-informative conditions. From this baseline, feature values are progressively introduced to quantify their contributions to model decision. To ensure physically meaningful inputs, we defined the baseline as the median vector computed across negative samples (non-STHP events, label = 0) in the training set. Specifically, for each of the 71 physical variables, we calculated the median value from all negative samples to construct a 71-dimensional reference input. This choice is motivated by two key considerations: (1) Unlike zero vectors or random noise—which often correspond to physically unrealistic atmospheric states (e.g., zero humidity or zero geopotential height)—the median of non-STHP samples represents a climatologically stable and meteorologically plausible reference regime. (2) By anchoring the baseline to observed non-convective conditions, the IG integration path traverses realistic atmospheric transitions (e.g., from stable to unstable stratification), thereby yielding physically interpretable attributions that reflect how deviations from typical “no-precipitation” environments contribute to STHP occurrence. This baseline represents typical non-convective environmental conditions, ensuring the integration path remains within meteorologically stable regimes. Compared to model outputs F ( x ′ ) derived from this baseline typically exhibit low probabilities, this approach aligns with the physical narrative of transitioning from “non-convective” to “potentially convective” conditions.

We implemented the IG method using the Python captum library,² which handles gradient computation and path integration through the Integrated Gradients function. For global feature importance assessment, we calculated absolute attribution values for all test samples and averaged them across the dataset. Results were visualized using swarm plots and bar charts to display attribution distributions and importance scores, while partial dependence plots were employed to analyze how interactions between key physical variables influence STHP identification.

All algorithms were evaluated on the same dataset randomly partitioned into training and test sets using an 8:2 ratio. The training set served for hyperparameter optimization and model development, while the test set was reserved exclusively for performance evaluation. To ensure reproducibility and minimize randomness in data partitioning, all experiments employed a fixed random seed (seed = 42). All models were trained and evaluated on identical data splits, ensuring fair and consistent comparative analysis. This study formulates STHP identification as a binary classification task, where samples with STHP occurrence are labeled positive and non-STHP samples as negative. Model performance was assessed using a confusion matrix (Table 3) and threat score (TS), supplemented by derived metrics including true positives (TP), false positives (FP), true negatives (TN), and false negatives (FN).

The TS and Matthews correlation coefficient (MCC) were calculated as primary evaluation metrics. The TS (Equation 8), widely adopted in the verification of meteorological event identification (Zhao, 2022) and assesses the accuracy of event-based identification by accounting for TP, FP, and FN. The MCC (Equation 9) provides a balanced measure of classification quality across both classes, representing the correlation between identified and observed binary outcomes (Chicco and Jurman, 2023; Stoica and Babu, 2024). The computational formulas are as follows:

The TS ranges from 0 to 1, where a value of 1 indicates perfect identification and 0 denotes complete identification failure. The MCC ranges from −1 to 1, where 1 represents perfect identification, 0 equates to random guessing, and −1 indicates total disagreement between identified and observed outcomes. Both metrics place particular emphasis on the accurate identification of STHP events, making them suitable for evaluating the rate of missed detections and false alarms of the models.

This study evaluates the performance of the proposed MSDNN for STHP identification against three benchmark models—LGB, RF, and DNN. Comprehensive performance metrics across both training and test datasets are systematically compared in Table 4.

On the training set, the MSDNN demonstrated exceptional capability, achieving the highest TS of 0.952 and MCC of 0.953, surpassing all comparative models. LGB (TS = 0.907, MCC = 0.906) and RF (TS = 0.755, MCC = 0.723) showed substantial performance gaps relative to the MSDNN. Confusion matrix analysis (Figure 4) further revealed that the MSDNN correctly identified 98.9% of positive samples (STHP events), exceeding the DNN (98.3%), LGB (97.5%), and RF (92.1%). For negative samples, the MSDNN attained a recognition rate of 96.5%, outperforming the DNN (94.5%), LightGBM (93.2%), and RF (80.0%). On the test set, although all models experienced performance degradation, the MSDNN maintained superior ability with a TS of 0.851 and MCC of 0.844. These results slightly exceeded those of the DNN (TS = 0.842, MCC = 0.832) and substantially outperformed LGB (TS = 0.808, MCC = 0.791) and RF (TS = 0.711, MCC = 0.663), indicating limited generalization capacity of the traditional ensemble methods in operational identification scenarios. Confusion matrix analysis (Figure 4) for the test set confirmed that the MSDNN correctly identified 95.3% of positive samples, marginally higher than the DNN (95.1%) and LGB (93.3%), and significantly exceeding RF (89.0%). For negative samples, the MSDNN achieved a recognition rate of 89.2%, again surpassing the DNN (88.2%), LGB (85.9%), and RF (77.1%).

The MSDNN consistently achieved optimal performance across both primary evaluation metrics (TS and MCC), demonstrating statistically significant advantages over traditional ensemble methods and a measurable improvement over the standard DNN architecture. Notably, the model exhibited enhanced stability and superior generalization capability in accurately identification both positive and negative cases, underscoring its practical potential for operational STHP identification.

Accurate identification of STHP requires models capable of discerning complex correlations among interconnected physical variables. This study proposes a MSDNN that demonstrates substantially superior performance compared to traditional machine learning models and conventional DNN. The MSDNN’s enhanced identification capability stems from its innovative architecture featuring five dedicated input branches for processing distinct meteorological variable categories. This design enables comprehensive feature extraction while preserving the unique physical properties of each variable type, with subsequent feature fusion creating highly discriminative representations. The “divide-and-conquer” approach effectively captures complex interdependencies and nonlinear interactions among diverse meteorological factors, significantly improving STHP identification.

Compared to traditional ensemble learning methods such as LGB and RF, the MSDNN demonstrates distinct advantages in handling STHP’s complex patterns. While ensemble methods can provide stable identification performance in certain scenarios, they often struggle to achieve sufficient accuracy when confronted with multi-source data and strong nonlinear relationships (Zhang et al., 2025). Existing research has established that DL architectures excel through hierarchical nonlinear transformations that automatically extract multi-level feature representations, providing inherent advantages in capturing complex meteorological processes (Meenal et al., 2021; Apaydın et al., 2022). In this study, although the benchmark DNN slightly underperformed relative to the MSDNN, it still surpassed both LGB and RF in identification capability, maintaining superior accuracy even under constrained computational resources. The MSDNN’s effectiveness is further supported by evidence from agricultural remote sensing, where similar multi-stream architectures have demonstrated robust multi-source data processing capabilities. For instance, Khan et al. (2025) achieved an R² of 0.86 for corn yield prediction using a two-stream DNN, significantly outperforming traditional methods and single-stream DNNs. Similarly, Li et al. (2025) developed a multi-source deep feature network for wheat chlorophyll estimation, attaining an R² of 0.85 with a 6.4% improvement over conventional DNN. These findings collectively confirm that multi-stream architectures can effectively capture inherent characteristics and interactive relationships among variables with distinct physical meanings, strongly validating our study conclusions.

This study applies the IG method to perform attribution analysis on the MSDNN model. The primary contribution lies not only in validating the discriminative logic of the model, but more importantly, in revealing complex nonlinear interactions and triggering mechanisms among environmental physical variables during STHP evolution. By combining model-derived feature attributions with classical convection theories, we elucidate the physical processes governing STHP occurrence.

The 700 hPa relative humidity (RH_700) emerges as the dominant moisture discriminator, displaying a marked nonlinear threshold effect. IG analysis indicates that intermediate RH_700 values (40–60%) strongly suppress STHP identification. In the Yangtze–Huaihe region, low mid-level humidity typically signals dry intrusion with northerly upper-level winds. When coupled with warm, moist low-level flows, this generates strong vertical wind shear that organizes convection but accelerates cell propagation, shortens precipitation duration, and restricts accumulation for extreme events. Dry intrusion further promotes evaporative cooling and entrainment, inhibiting upscale convective growth. Conversely, RH_700 exceeding 80% sharply increases positive attribution: a deep moist layer supplies abundant vapor, reduces mid-level dissipation, preserves warm-core structures, and supports transition to efficient heavy precipitation. These results align with Kang and Ebtehaj (2025), who found moisture variables dominant in convective discrimination, and Yousefnia et al. (2025), who emphasized low-level RH for thunderstorm initiation. The specific humidity threshold of 9 g/kg similarly marks a high-contribution regime, reflecting Clausius–Clapeyron constraints requiring warm temperatures to sustain such moisture. The model’s capture of this “high temperature–high humidity” coupling matches observational thresholds (RH_700 of 76%; specific humidity of 8.36 g/kg) from Shen et al. (2016).

WVV contribute 45% of total attribution, highlighting the pivotal role of large-scale synoptic forcing in regional STHP. 500 hPa vertical speed is the key factor, with negative values yielding strong positive contributions. This reflects mesoscale control by favorable large-scale circulation. Quasi-geostrophic theory positions the 500 hPa level near nondivergence, where ascent is weak under balance; only deep troughs and positive vorticity advection induce organized uplift. These findings corroborate Retsch et al. (2022) on vertical velocity in deep convection and Fan et al. (2024) SHAP-based sensitivities. Notably, the model integrates southwesterly winds across levels (950–200 hPa), capturing low-level jet convergence and upper-level divergence as a “dynamic pump” that compensates for modest CAPE and enables extreme precipitation under weakly unstable conditions. This contrasts with thermodynamically focused studies (e.g., Zhang et al., 2019) and demonstrates MSDNN’s capacity to autonomously extract three-dimensional dynamical structures.

Attributions related to thermodynamic conditions and vertical structure reveal a delicate balance between energy release and precipitation efficiency. CAPE contributions peak near 3,000 J/kg before declining, suggesting microphysical limits: excessive CAPE may drive intense updrafts favoring ice-phase loading and reducing warm-rain efficiency. A high 700 hPa pseudo-equivalent potential temperature ( θ se ≈ 340 K) inhibits by flattening the vertical θ se gradient and suppressing low-level instability, consistent with eastern China observations (Qiu and He, 2013). Small vertical separation (<500 m) between the lifting condensation level (LCL) and convective condensation level (CCL) promotes STHP by enabling efficient warm-rain processes through high near-surface moisture and low trigger heights.

In summary, IG analysis not only confirms MSDNN’s reliability in STHP identification but also quantitatively reconstructs a physically coherent convective scenario aligned with the “ingredients-based” framework (Doswell et al., 1996). Deep moist layers with warm backgrounds mitigate entrainment (moisture constraint); large-scale troughs and dynamic pumping disrupt balance (dynamic constraint); and moderate energy with optimized stratification (thermodynamic constraint) collectively define the conditions for extreme precipitation.

This study presents an innovative framework that integrates a MSDNN with IG interpretability analysis, establishing a novel methodology for STHP identification that effectively balances identification performance with physical interpretability. The principal strengths of this research are manifested in two key aspects. First, the MSDNN architecture is explicitly designed to align with established physical principles of severe convective weather. By categorizing 71 environmental physical variables into five specialized input streams—moisture conditions, dynamic parameters, thermal instability, comprehensive indices, and vertical level features—the model facilitates hierarchical extraction of features from heterogeneous meteorological elements. This design enables dedicated sub-networks to capture intrinsic patterns within specific physical processes, while a subsequent integration layer models cross-process nonlinear interactions. Validation result confirms that the MSDNN surpasses both conventional ensemble learning methods (LGB, RF) and single-stream DNN, demonstrating that this physics-informed partitioning strategy more effectively utilizes multi-source meteorological information and enhances discrimination capability for complex convective systems. Second, the successful application of the integrated gradients method provides an effective solution to the “black-box” problem of DL-based meteorological models. Beyond quantifying the marginal contributions of individual physical variables to model outputs, IG analysis reveals nonlinear threshold effects of key factors. These findings offer data-driven new perspectives for understanding the initiation mechanisms of STHP. Through feature attribution analysis, we establish a quantitative link between model identification and physical processes, enhancing the credibility of the model in operational applications.

However, several limitations warrant consideration. First, the station-based 1D feature representation, while enabling direct physical interpretation, does not encode spatial neighborhood dependencies critical for capturing mesoscale convective organization and frontal structures—capabilities inherent to 2D convolutional architectures like U-Net. Second, the ERA5 reanalysis data employed in this study, while comprehensive, lacks the spatial resolution necessary to fully resolve fine-scale structures of mesoscale and microscale convective systems. This scale discrepancy consequently leads to elevated false alarm rates in identifying isolated severe convection events, particularly those initiated by subgrid-scale processes such as boundary layer convergence lines and localized frontal zones. Furthermore, the current feature representation system inadequately captures topographic forcing effects, despite the demonstrated significance of Anhui Province’s complex terrain—including the Dabie Mountains and southern mountainous regions—in modulating precipitation distribution through mechanical lifting and moisture channeling mechanisms. The existing environmental physical variables prove insufficient in characterizing these topographically driven dynamic processes. Finally, the current framework utilizes reanalysis-based diagnostics to identify STHP precursors, which provides a retrospective understanding of critical atmospheric conditions but does not directly evaluate predictive skill at operational forecast lead times. Future research will address these limitations through three primary directions: (1) We will develop grid-based multi-branch architectures that combine spatial modeling techniques—such as U-Net for local feature extraction and Transformer for long-range dependency capture—with existing physics-based variable streams, thereby integrating spatial neighborhood correlations into the identification framework. (2) We will prioritize the integration and evaluation of high-resolution regional reanalysis products to better resolve convective-scale atmospheric structures. Concurrently, we plan to incorporate a dedicated topographic input branch within the MSDNN framework, integrating parameters such as elevation, slope gradient, aspect, and underlying surface characteristics. (3) To facilitate operational deployment, a transfer learning strategy (Ham et al., 2019) will be employed. The MSDNN model pre-trained on ERA5 reanalysis data will serve as the initialization, with fine-tuning conducted using bias-corrected operational NWP fields. These enhancements are expected to strengthen the coupling between MSDNN and numerical weather prediction systems, advancing severe convective weather forecasting toward greater precision and enhanced physical interpretability.