A method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling

Abstract

The application discloses a kind of based on multi-factor coupling's hill area extreme rainstorm event prediction and analysis method, comprising: step S1 data collection and pre-processing, the data after pre-processing is divided into training dataset and test dataset, step S2 multi-factor feature extraction, extract main meteorological characteristics, calculate topographic index to represent the influence of topography on precipitation, step S3 key influence factor screening, step S4 prediction model construction and training, adopt support vector machine algorithm to build extreme rainstorm event prediction model, step S5 prediction result evaluation, including accuracy, recall rate, precision, balance the comprehensive index F1 value of the recall rate and precision.The method of the application reduces the redundancy and noise of data, improves the efficiency and accuracy of the prediction model, can more accurately capture the complex relationship between multi-factors, thereby improving the prediction accuracy of extreme rainstorm event.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of meteorological disaster prediction, and particularly relates to a mountainous area extreme rainstorm event prediction and analysis method based on multi-factor coupling. BACKGROUND

[0002] Under the background of global climate change, extreme rainstorm events occur frequently, causing serious natural disasters such as landslides, debris flows, and floods in mountainous areas, which pose a great threat to people's life and property safety and the ecological environment.

[0003] AI and machine learning are widely used in weather prediction, and existing research has extensively explored how machine learning can be applied to rainfall prediction and disaster prevention and mitigation. However, existing rainstorm prediction methods mostly focus on the analysis of a single meteorological factor or a few factors, especially in mountainous and complex terrain conditions. The existing technology fails to fully consider the coupling effect between the complex topography, atmospheric circulation, and water vapor transport in mountainous areas, resulting in insufficient accuracy and reliability in predicting extreme rainstorm events.

[0004] Therefore, it is of great practical significance to develop a mountainous area extreme rainstorm event prediction and analysis method that can comprehensively consider the influence of multiple factors. SUMMARY

[0005] To solve the above problems, the present application proposes a mountainous area extreme rainstorm event prediction and analysis method based on multi-factor coupling, which uses a combination of principal component analysis, random forest, and support vector machine to build a prediction model by comprehensively analyzing multiple meteorological and geographical factors, to improve the prediction accuracy of extreme rainstorm events in mountainous areas.

[0006] The present application is implemented as follows: a mountainous area extreme rainstorm event prediction and analysis method based on multi-factor coupling, comprising the following steps:

[0007] Step S1, data collection and preprocessing: collect ERA5 dataset of the European Centre for Medium-Range Weather Forecasts, measured data of meteorological stations in the study area, and DEM digital elevation model data reflecting the topographic features of the study area; preprocess the collected data, including quality control and standardization processing, the quality control includes using a hybrid strategy of spatiotemporal coupling dynamic threshold method + multivariate joint test to remove outliers and fill missing values, and dividing the preprocessed data into training dataset and test dataset;

[0008] Step S2, multi-factor feature extraction: use principal component analysis method for dimensionality reduction, extract main meteorological features, and calculate terrain index to represent the influence of terrain on precipitation;

[0009] Step S3, Key Influencing Factor Screening: Use the random forest algorithm to rank the importance of the extracted multi-factor features and select the features with higher importance scores as the key influencing factors of extreme rainstorm events;

[0010] Step S4, Prediction Model Construction and Training: A prediction model for extreme rainstorm events is constructed using the support vector machine algorithm; the key influencing factors of the selected extreme rainstorm events are used as the input features of the prediction model, and "whether the extreme rainstorm event has occurred" is used as the output label of the prediction model. The prediction model is trained until the model training is completed.

[0011] Step S5, Evaluation of prediction results: The extreme rainstorm event prediction model trained in step S4 is evaluated using the test dataset described in step S1. The evaluation metrics include accuracy, recall, precision, and the F1 score, which is a comprehensive metric balancing the recall and precision.

[0012] Furthermore, the specific steps for removing outliers in step S1 include:

[0013] The first step, univariate spatiotemporal dynamic threshold detection: construct a spatiotemporal sliding window for each meteorological variable and calculate the dynamic threshold Q. upper (t, s):

[0014]

[0015] In the formula, Q 75 (t, s) represents the 75th percentile within the spatiotemporal window, IQR(t, s) represents the interquartile range, and IQR(t, s) = Q 75 (t, s)-Q 25 (t, s), Q 25 (t, s) represents the 25th percentile within the spatiotemporal window; Δh represents the mean elevation difference between the weather station and its neighboring area, h ref For reference elevation parameters, k(h) is the topographic adjustment factor, h is the elevation of the meteorological station, and h meam σ represents the average elevation of the study area. h Standard deviation of elevation;

[0016] For suspected outliers exceeding the threshold, the time series curvature is used for analysis. Verification, the formula is:

[0017]

[0018] A value is considered an outlier if it meets the following conditions:

[0019] and

[0020] In the formula, t is the time step at the current moment, and xt Let x be the observed value of the meteorological variable at time t. t+1 x is the observation at the next time step after time t. t-1 The observation value is the one time step prior to time t. For the standard deviation of curvature, This is the seasonal average.

[0021] The second step is multivariate joint anomaly detection: calculating the multivariate Mahalanobis distance D. M :

[0022]

[0023] In the formula, X is the multivariate observation vector, μ is the mean vector, T is the transpose operator of a matrix or vector, and Σ is the covariance matrix;

[0024] Introducing terrain weight correction threshold D threshold :

[0025]

[0026] In the formula, V represents the number of variables. TI is the 99th quantile of a chi-square distribution with V degrees of freedom, and TI is the topographic index of the current location. max The maximum topographic index of the study area is given by 'a', where 'a' is the catchment area and 'b' is the slope.

[0027] If D is satisfied M >D threshold This is then identified as a multivariate coupling outlier;

[0028] Step 3: Spatiotemporal density clustering verification: Calculate the composite spatiotemporal distance metric d st :

[0029]

[0030] In the formula, d space Let d be the spatial Euclidean distance. time d is the time interval. max To determine the maximum Euclidean distance between all pairs of spatial points within the study region, ΔT max To study the maximum time span within the time range, Corr(X, Y) is the Pearson correlation coefficient of the variable trend, and α, β, γ are adjustable weight coefficients, where α+β+γ=1;

[0031] Calculate the adaptive density parameters:

[0032] ∈=μ d +2σ d

[0033]

[0034] In the formula, ∈ represents the neighborhood radius, and μ d Let σ be the mean of the distance distribution. d N is the distance standard deviation, MinPts is the minimum number of cluster points, and N is the minimum number of cluster points. total This represents the total number of samples.

[0035] If a data sample point does not belong to any density cluster and is not covered by more than 80% of its neighboring points, it is considered an outlier.

[0036] The fourth step is physical mechanism inversion verification: Water vapor flux is used to verify the outliers to be removed. Verification, the formula is:

[0037]

[0038] If the following conditions are met, the event is considered a real extreme event, and the data will be retained:

[0039] (Historical extreme values ​​of the study area)

[0040] In the formula, Q 95 The value represents the 95th percentile of the historical extreme values ​​for the study area, where q represents specific humidity, i.e., the mass of water vapor per unit mass of moist air. This is a three-dimensional wind speed vector, where x represents the east-west direction, y represents the north-south direction, z represents the vertical direction, u represents the wind speed in the east-west direction (east is positive), v represents the wind speed in the north-south direction (north is positive), and w represents the wind speed in the vertical direction (upward is positive).

[0041] Furthermore, the filling of missing values ​​in step S1 specifically involves using linear interpolation to fill in the missing values, assuming the original data sequence is x1, x2, ..., x... m The linear interpolation formula is:

[0042]

[0043] In the formula, x i For missing values, x i-k For x i Previous adjacent valid data, x i+l For x i The next valid data points are k and l, where k and l represent the distance from the missing value position i to the previous valid data point x. i-k The next adjacent valid data x i+l The number of intervals.

[0044] Furthermore, the standardization process described in step S1 specifically involves transforming the values ​​of each variable into a standard normal distribution with a mean of 0 and a standard deviation of 1, in order to eliminate the influence of different variable units. The standardization formula is as follows:

[0045]

[0046] In the formula, is the value after standardization of the j-th variable value of the i-th sample, and x ij is the j-th variable value of the i-th sample, is the mean of the j-th variable, and s j is the standard deviation of the j-th variable.

[0047] Furthermore, the multi-factor feature extraction described in step S2 specifically includes the following sub-steps:

[0048] Step S21, calculate the covariance matrix S: Let the measured meteorological data of the meteorological stations in the training data set described in step S1 be an n×p data matrix X * , where n is the number of samples and p is the total number of meteorological features; calculate the X * covariance matrix S:

[0049]

[0050] In the formula, T is the transpose operator of a matrix or vector, and X *T is the transpose matrix of X * , and the covariance matrix S is a p×p real symmetric matrix;

[0051] Step S22, eigenvalue decomposition: For the covariance matrix S obtained in step S21, there exist an orthogonal matrix E and a diagonal matrix Λ such that:

[0052] S = EΛE T

[0053] where the diagonal elements of the diagonal matrix Λ are the eigenvalues λ1 ≥ λ2 ≥ … ≥ λ p of the covariance matrix S, and the column vectors of the orthogonal matrix E are the corresponding unit orthogonal eigenvectors e1, e2, …, e p ;

[0054] Step S23, dimensionality reduction processing: Select the eigenvectors e1, e2, …, e k corresponding to the first k eigenvalues as column vectors respectively to form a projection matrix P:

[0055] P = (e1, e2, …, e k )

[0056] where the value range of k is 1 ≤ k < p; when selecting the first k eigenvalues, ensure that the cumulative variance contribution rate is greater than or equal to the first preset threshold, and the formula for the cumulative variance contribution rate Cumulative Variance Ratio is:

[0057]

[0058] In the formula, λ i is the i-th eigenvalue sorted from large to small;

[0059] The matrix X of the measured meteorological station data in the training data set in step S1 is projected into a low-dimensional space through the projection matrix P * to obtain the principal component matrix Y:

[0060] Y = XP * P

[0061] Finally, the first k main meteorological features are extracted;

[0062] Step S24, calculate the topographic index: Calculate the topographic index TI for the digital elevation model data of the training data set in step S1:

[0063]

[0064] In the formula, a is the catchment area and b is the slope.

[0065] Further, the screening of the key influencing factors in step S3 specifically includes the following sub-steps:

[0066] Step S31, calculate the overall importance: First, use a random forest to construct multiple decision trees to classify or regress the data. Each decision tree randomly selects a part of the samples and features for training. After training, evaluate its overall importance by calculating the average impurity reduction of each feature in all decision trees. Let the importance of feature j in the t-th decision tree be I j，t , and the formula for the overall importance of feature j is:

[0067]

[0068] In the formula, I j is the overall importance of feature j, α is the number of decision trees, and I j，t is the importance of feature j in the t-th decision tree;

[0069] Step S32, select the key influencing factors: Sort the multi-factor features extracted in step S2 according to the overall importance from large to small, and select the top m features with higher overall importance as the key influencing factors of extreme rainstorm events; where the value range of m is 1 ≤ m < q, and q is the total number of multi-factor features extracted in step S2; when selecting the top m features with higher overall importance, ensure that the cumulative importance is greater than or equal to the second preset threshold. The formula for the cumulative importance Cumulative Importance is:

[0070]

[0071] In the formula, I j The overall importance of feature j, sorted from largest to smallest.

[0072] Furthermore, the support vector machine algorithm described in step S4 is as follows: For the linearly separable case, the formula for solving the optimization problem using support vector machine is:

[0073]

[0074] sty i (w T φ(x i )+b)≥1, i=1, 2,…,n

[0075] In the formula, w is the weight vector, b is the bias term, and x i For the input sample, φ(x) i Let y be a function that maps input samples to a high-dimensional feature space. i y represents the class label of the sample. i ∈{-1, 1};

[0076] For the linearly inseparable case, the formula for solving the optimization problem using support vector machines is:

[0077]

[0078] sty i (w T φ(x i )+b)≥1-ξ i ξ i ≥0, i=1,2,…,n

[0079] In the formula, ξ i C is a slack variable, and C is a penalty factor.

[0080] Furthermore, the specific process of training the prediction model until the model training is completed in step S4 is as follows: during the training process, grid search is used to traverse the combination of kernel function parameters, and cross-validation is used to optimize the model parameters, and the optimal penalty factor C and kernel function parameters are selected.

[0081] The completion of model training specifically includes the following criteria:

[0082] The optimization algorithm meets the convergence conditions: the change in the objective function between adjacent iterations is less than the third preset threshold; the training samples satisfy the Karush-Kuhn-Tucker conditions.

[0083] Preset upper limit for the number of iterations: If the optimization algorithm fails to converge within a fixed number of iterations, it will be forcibly stopped to avoid infinite loops;

[0084] Stable performance on the test dataset: If the model performance no longer improves after a certain number of consecutive rounds on the cross-validation or the test dataset, training will be terminated early to prevent overfitting.

[0085] Parameter optimization complete: Training is considered complete when the optimal penalty factor C and kernel function parameters are found through grid search, and the model's solution on the training set satisfies the constraints of the optimization problem.

[0086] Furthermore, the formula for calculating the accuracy in step S5 is as follows:

[0087]

[0088] In the formula, TP is the number of true positives, TN is the number of true negatives, FP is the number of false positives, and FN is the number of false negatives; the range of Accuracy is [0, 1], where 0 indicates that the prediction of the extreme rainstorm event is completely wrong, and 1 indicates that the prediction of the extreme rainstorm event is completely correct.

[0089] Furthermore, the F1 value mentioned in step S5 is:

[0090]

[0091] In the formula, the F1 value represents the harmonic mean of the two indicators, Precision and Recall, indicating that the two indicators are equally important; the range of the F1 value is [0, 1], where 0 means that at least one of Precision and Recall is 0, and 1 means that both Precision and Recall are 1.

[0092] The beneficial effects of this invention are as follows: It comprehensively considers the impact of various meteorological and geographical factors on extreme rainstorm events in hilly areas. Through multi-factor feature extraction and key influencing factor screening, it reduces data redundancy and noise, improving the efficiency and accuracy of the prediction model. Employing advanced data analysis and machine learning algorithms such as principal component analysis, random forest, and support vector machine, it can more accurately capture the complex relationships between multiple factors, thereby improving the prediction accuracy of extreme rainstorm events. Rigorous data preprocessing and model evaluation ensure the reliability and stability of the method, providing strong technical support for meteorological disaster early warning and disaster prevention and mitigation work in hilly areas.

[0093] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. Attached Figure Description

[0094] Figure 1 This is a flowchart of the method steps of the present invention. Detailed Implementation

[0095] Example 1:

[0096] This embodiment provides a method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling, such as... Figure 1 As shown, it includes the following steps:

[0097] Step S1, Data Collection and Preprocessing:

[0098] We collected the ERA5 dataset from the European Centre for Medium-Range Weather Forecasts (ECMWF), including hourly data on pressure levels from 1940 to present and monthly averaged data from 1950 to present for ERA5-Land. We also acquired measured data from meteorological stations in the study area, including meteorological elements such as precipitation, temperature, air pressure, wind speed, and wind direction. In addition, we collected digital elevation model (DEM) data reflecting the topographic features of the study area.

[0099] Use Python's pandas library to read data and preprocess the collected data, including quality control and standardization.

[0100] (1) Regarding data quality control, a hybrid strategy of spatiotemporal coupling dynamic thresholding and multivariate joint testing is adopted for outlier detection. The specific steps include:

[0101] Step 1, Univariate Spatiotemporal Dynamic Threshold Detection:

[0102] Using the terrain-adaptive sliding quantile method, a spatiotemporal sliding window is constructed for each meteorological variable (such as precipitation), and the dynamic threshold is calculated using the following formula:

[0103]

[0104] In the formula, Q 75 (t, s) represents the 75th percentile within the spatiotemporal window, IQR(t, s) represents the interquartile range, and IQR(t, s) = Q 75 (t, s)-Q 25 (t, s), Q 25 (t, s) represents the 25th percentile within the spatiotemporal window; Δh represents the mean elevation difference between the weather station and its neighboring area, h re′ For reference elevation parameters, k(h) is the topographic adjustment factor, h is the elevation of the meteorological station, and h meam σ represents the average elevation of the study area. h This represents the standard deviation of elevation.

[0105] Using the second derivative of the mutation point for verification, suspected outliers exceeding the threshold are analyzed by time series curvature. Verification, the formula is:

[0106]

[0107] A value is considered an outlier if it meets the following conditions:

[0108] and

[0109] In the formula, t is the time step or time index of the current moment (e.g., if the data is observed hourly, t corresponds to the t-th hour), x t Let x be the observed value of a meteorological variable at time t (e.g., precipitation). t+1 x represents the observation at the next time step after time t (e.g., precipitation at hour t+1). t-1 This refers to the observation value at the previous time step before time t (e.g., the precipitation at hour t-1). For the standard deviation of curvature, This is the seasonal average.

[0110] The second step is multivariate joint anomaly detection (terrain-corrected Mahalanobis distance):

[0111] Calculate the multivariable Mahalanobis distance D M :

[0112]

[0113] In the formula, X is the multivariate observation vector, μ is the mean vector, T is the transpose operator of a matrix or vector (i.e., the transpose symbol in linear algebra), and Σ is the covariance matrix.

[0114] Calculate the terrain weight correction threshold D threshold :

[0115]

[0116] In the formula, V represents the number of variables. TI is the 99th quantile of a chi-square distribution with V degrees of freedom, and TI is the topographic index of the current location. max The maximum topographic index of the study area is given by denoted as a, where a is the catchment area and b is the slope.

[0117] If D is satisfied M >D threshold It is then determined to be a multivariate coupling outlier.

[0118] The third step is the verification of spatiotemporal density clustering (an improved version of the ST-DBSCAN algorithm):

[0119] Calculate the composite spatiotemporal distance metric dst :

[0120]

[0121] In the formula, d space Let d be the spatial Euclidean distance. time d is the time interval. max To determine the maximum Euclidean distance between all pairs of spatial points within the study region, ΔT max To study the maximum time span within the time range, Corr(X, Y) is the Pearson correlation coefficient for the trend of the variables, and α, β, and γ are adjustable weighting coefficients, with α + β + γ = 1. In this embodiment, α = 0.4, β = 0.3, and γ = 0.3 are used.

[0122] Calculate the adaptive density parameters:

[0123] ∈=μ d +2σ d

[0124]

[0125] In the formula, ∈ represents the neighborhood radius, and μ d Let σ be the mean of the distance distribution. d N is the distance standard deviation, MinPts is the minimum number of cluster points, and N is the minimum number of cluster points. total This represents the total number of samples.

[0126] If a data sample point is an isolated point, that is, it does not belong to any density cluster and is not covered by more than 80% of its neighboring points, it is considered an outlier.

[0127] Step 4: Physical mechanism inversion verification:

[0128] Water vapor flux was used to measure the outliers to be removed. Verification, the formula is:

[0129]

[0130] If the following conditions are met, the event is considered a real extreme event, and the data will be retained:

[0131] (Historical extreme values ​​of the study area)

[0132] In the formula, Q 95 The value represents the 95th percentile of the historical extreme values ​​for the study area. q represents specific humidity, which is the mass of water vapor per unit mass of moist air, characterizing air humidity. The three-dimensional wind speed vector consists of three directional components: u, v, and w. Let x be the east-west direction (horizontal direction), y be the north-south direction (horizontal direction), and z be the vertical direction (height or altitude). Then u is the wind speed in the east-west direction (east is positive), v is the wind speed in the north-south direction (north is positive), and w is the wind speed in the vertical direction (upward is positive).

[0133] (2) After an outlier is detected, if the proportion of outliers is small (the proportion of outliers in this embodiment is <5%) and does not affect the overall distribution and the amount of data is sufficient, then these outliers are removed.

[0134] (3) For missing values, linear interpolation is used for filling. The specific method is as follows:

[0135] Let the original data sequence be x1, x2, ..., x m The linear interpolation formula is:

[0136]

[0137] In the formula, x i For missing values, x i-k For x i Previous adjacent valid data, x i+l For x i The next valid data points are k and l, where k and l represent the distance from the missing value position i to the previous valid data point x. i-k The next adjacent valid data x i+l The number of intervals, i.e., it is necessary to skip k positions forward and l positions backward to find the nearest valid data.

[0138] (4) Because subsequent steps require calculating the covariance matrix, which is sensitive to the dimensions of variables, if the dimensions of the original data features differ greatly, such as temperature (°C or °F) and rainfall (mm), directly calculating the covariance matrix will amplify the impact of high-dimensional features on the results. Therefore, it is necessary to standardize the data first. The specific method for data standardization is to use the StandardScaler class in the scikit-learn library to transform the values ​​of each variable into a standard normal distribution with a mean of 0 and a standard deviation of 1, so as to eliminate the influence of different variable dimensions (Z-score standardization).

[0139] The standardization formula is:

[0140]

[0141] In the formula, Let x be the standardized value of the j-th variable in the i-th sample. ij Let j be the value of the j-th variable in the i-th sample. Let s be the mean of the j-th variable.j Let be the standard deviation of the j-th variable. After standardization, the covariance matrix is ​​equivalent to the correlation coefficient matrix of the original data.

[0142] (5) After preprocessing, the dataset is divided into training dataset and test dataset according to a certain ratio.

[0143] Step S2, Multi-factor Feature Extraction:

[0144] The measured meteorological data from the weather stations in the training dataset obtained in step S1 are used to perform dimensionality reduction processing using the PCA class (Principal Component Analysis) method from the scikit-learn library to extract the main meteorological features. This specifically includes the following sub-steps:

[0145] Step S21, calculate the covariance matrix S: Let the measured meteorological station data in the training dataset of step S1 be an n×p data matrix X. * Where n is the number of samples and p is the total number of meteorological features; the calculated X * Covariance matrix S:

[0146]

[0147] In the formula, T is the transpose operator for a matrix or vector (i.e., the transpose symbol in linear algebra), and X... *T For X * The transpose of the matrix, and the covariance matrix S is a p×p real symmetric matrix (S T =S).

[0148] Step S22, Eigenvalue decomposition: According to linear algebra theory, a real symmetric matrix S can be orthogonally diagonalized, that is, there exists an orthogonal matrix E(E T E = I) and a diagonal matrix Λ, such that:

[0149] S = EΛE T

[0150] The diagonal elements of the diagonal matrix Λ are the eigenvalues ​​of the covariance matrix S: λ1≥λ2≥…≥λ p The column vectors of the orthogonal matrix E are the corresponding unit orthogonal eigenvectors e1, e2, ..., e p Numerical methods (such as the Jacobi method and the QR algorithm) are used to solve for eigenvalues ​​and eigenvectors.

[0151] It is important to note that NumPy returns eigenvalues ​​of real symmetric matrices in ascending order, and this order needs to be reversed to satisfy λ1≥λ2≥…≥λ p The eigenvectors of a real symmetric matrix are naturally orthogonal. And it is a unit vector (‖e i‖ = 1), which can be verified by calculating the dot product through NumPy.

[0152] Step S23, dimensionality reduction processing: Select the eigenvectors e1, e2, …, e corresponding to the first k eigenvalues k as column vectors respectively to form the projection matrix P:

[0153] P = (e1, e2, …, e k )

[0154] where the value range of k is 1 ≤) < p. Since at least 1 principal component needs to be retained and it must be less than the original feature number p to achieve the purpose of dimensionality reduction (if k = p, there is no dimensionality reduction).

[0155] When selecting the first k eigenvalues, ensure that the cumulative variance contribution rate is greater than or equal to the first preset threshold (85% is used in this embodiment, and it can also be changed to commonly used thresholds in the industry such as 90%, 95%, etc.). The formula for the cumulative variance contribution rate Cumulative Variance Ratio is:

[0156]

[0157] In the formula, λ i is the i-th eigenvalue sorted from large to small.

[0158] When the cumulative variance contribution rate reaches the preset 85%, it is considered that the first k principal components have retained the main information of the original data and can be used as the "main meteorological features". In traditional methods, the principal components with eigenvalues ≥ 1 are often selected because the components with eigenvalues less than 1 carry less information than a single original feature. However, in practical applications, if the cumulative variance contribution rate of the principal components with eigenvalues ≥ 1 does not reach the threshold, additional subsequent principal components need to be supplemented until the threshold requirement is met. If the efficiency of the model is pursued, the threshold can be appropriately reduced (such as 80%) to reduce the value of k; if more detailed information needs to be retained, the threshold can be increased (such as 95%) to increase the value of k. Finally, the optimal k is determined through a scree plot or cross-validation to ensure that the data after dimensionality reduction can retain sufficient information and reduce redundancy.

[0159] Project the matrix X of the measured data of the meteorological station through the projection matrix P * onto the low-dimensional space to obtain the principal component matrix Y:

[0160] Y = X * P

[0161] Finally, the first k main meteorological features are extracted.

[0162] Step S24, Calculate the topographic index: For the digital elevation model data (topography) in the training dataset obtained in Step S1, use GIS software (such as ArcGIS) to calculate the topographic index to represent the impact of topography on precipitation. The formula is:

[0163]

[0164] In the formula, TI is the topographic index (TopographicIndex), a is the catchment area, and b is the slope.

[0165] Step S3, Screen key influencing factors:

[0166] Use the RandomForestClassifier class in the scikit-learn library to build a random forest model. Use the random forest algorithm to rank the importance of the multi-factor features extracted in Step S2. According to the importance ranking results, select the features with higher importance scores as the key influencing factors of extreme rainstorm events. Specifically, it includes the following sub-steps:

[0167] Step S31, Calculate the overall importance: First, use the random forest to build multiple decision trees to classify or regress the data. Each decision tree randomly selects a part of the samples and features for training. After training, evaluate its overall importance by calculating the average decrease in impurity (Gini index) of each feature in all decision trees.

[0168] Let the importance of feature j in the t-th decision tree be I j，t , and the formula for the overall importance of feature j is:

[0169]

[0170] In the formula, I j is the overall importance of feature j, α is the number of decision trees, and I j，t is the importance of feature j in the t-th decision tree.

[0171] Step S32, Select key influencing factors: Sort the multi-factor features extracted in Step S2 in descending order of overall importance, and select the top m features with higher importance scores as the key influencing factors of extreme rainstorm events.

[0172] Among them, the value range of m is 1 ≤ m < q, where q is the total number of multi-factor features extracted in Step S2; because at least 1 key influencing factor must be retained and it must be less than the total number of features q to achieve the screening purpose (if k = q, there is no screening).

[0173] When selecting the top m features with high overall importance, ensure that the cumulative importance is greater than or equal to a second preset threshold (e.g., 80%, 85%, or 90%, this embodiment uses 85%). The formula for cumulative importance is:

[0174]

[0175] In the formula, I j The overall importance of feature j, sorted from largest to smallest.

[0176] When the cumulative importance reaches a preset 85%, the top m features are considered to cover the key factors with the greatest impact on extreme rainstorm events and can be used as input features for subsequent model training. Alternatively, a minimum threshold for the importance score of a single feature can be set (e.g., greater than 5%), selecting only the top m features with importance scores above this threshold to avoid including features with weak contributions to the model. Furthermore, expertise in extreme rainstorms in hilly areas can be incorporated to ensure that key geographical features such as topographic indices are prioritized; even if their importance scores are slightly lower, the value of m can be adjusted manually.

[0177] Alternatively, the optimal m can be determined through cross-validation / grid search: train and test the model with different m values ​​(e.g., m = 5, 10, 15), and select the m that maximizes prediction accuracy. This balances model complexity with generalization ability, avoiding underfitting due to insufficient features or overfitting due to excessive features.

[0178] Step S4, Prediction Model Construction and Training:

[0179] (1) Use the SVC class in the scikit-learn library, i.e., the Support Vector Machine (SVM) algorithm, to construct a prediction model for extreme rainstorm events. The basic idea of ​​SVM is to find an optimal hyperplane in the feature space such that samples of different classes can be separated to the greatest extent. For the linearly separable case, the goal of SVM is to solve the following optimization problem:

[0180]

[0181] sty i (w T φ(x i ()+b)≥1, i=1,2,…,n

[0182] In the formula, w is the weight vector, b is the bias term, and x i For the input sample, φ(x) i Let y be a function that maps input samples to a high-dimensional feature space. i y represents the class label of the sample. i ∈{-1, 1}.

[0183] For the linearly inseparable case, a slack variable ξ is introduced. i With ≥0 and a penalty factor C, the optimization problem becomes:

[0184]

[0185] sty i (w T φ(x i )+b)≥1-ξ i ξ i ≥0, i=1,2,…,n

[0186] In the formula, ξ i Let C be a slack variable and C be a penalty factor.

[0187] (2) Using the key influencing factors of extreme rainstorm events selected in step S3 as input features of the prediction model, and "whether an extreme rainstorm event has occurred" as the output label of the prediction model, the prediction model is trained. The specific process is as follows:

[0188] During training, GridSearchCV is used to traverse the kernel function parameter combinations using grid search, and cross-validation (CV) is used to optimize the model parameters, selecting the optimal penalty factor C and kernel function parameters (such as gamma).

[0189] The completion of model training includes the following criteria:

[0190] SVM finds the optimal hyperplane by solving a convex optimization problem. The core criterion for successful training is that the optimization algorithm satisfies the convergence condition: the objective function between adjacent iterations... The change is less than the third preset threshold (10 is used in this embodiment). -6 The training samples satisfy the Karush-Kuhn-Tucker conditions, that is, the Lagrange multipliers of the support vectors are non-zero, and the constraints of the non-support vectors are strictly satisfied.

[0191] Preset maximum number of iterations: If the optimization algorithm fails to converge within a fixed number of iterations (e.g., 1000), it will be forcibly stopped to avoid infinite loops. This is set via the `max_iter` parameter in scikit-learn.

[0192] Stable performance on the test dataset: If the model performance no longer improves after a certain number of rounds (e.g., 10 rounds) on the cross-validation or test dataset, training is terminated early to prevent overfitting.

[0193] Parameter optimization complete: When the optimal penalty factor C and kernel function parameters are found through grid search, and the model's solution on the training set satisfies the constraints of the optimization problem (such as the slack variable constraint y in the case of linear inseparability). i (w T φ(x i )+b)≥1-ξ i Once completed, the training is considered finished.

[0194] Step S5, Evaluation of prediction results:

[0195] The extreme rainstorm event prediction model trained in step S4 was evaluated using the test dataset obtained in step S1. Evaluation metrics included accuracy, recall, precision, and F1 score to assess model performance. Accuracy reflects the overall correctness of predictions and the model's overall ability to distinguish between extreme and non-extreme rainstorm events. Recall reflects the model's ability to capture positive examples (i.e., extreme rainstorm events) to avoid missing predictions (leading to incomplete disaster warnings and high risk). Precision reflects the reliability of positive predictions to reduce false positives (potentially leading to over-preparation and increased social costs). The F1 score is a comprehensive metric balancing recall and precision, aiming to achieve a balance between missed and false positives. The evaluation metrics were calculated using the `classification_report` function from the scikit-learn library.

[0196] The formula for calculating accuracy is:

[0197]

[0198] In the formula, TP is the number of true positives (predicted as positive and actually positive), TN is the number of true negatives (predicted as negative and actually negative), FP is the number of false positives (predicted as positive but actually negative), and FN is the number of false negatives (predicted as negative but actually positive). The accuracy range is [0,1], where 0 indicates that the prediction of the extreme rainstorm event is completely wrong (all samples are misclassified), and 1 indicates that the prediction of the extreme rainstorm event is completely correct (no misclassification).

[0199] Recall rate is:

[0200]

[0201] In the formula, the range of Recall is [0,1]. 0 indicates that no positive examples were correctly identified (TP=0, all positive examples were misclassified as negative examples), and 1 indicates that all positive examples were correctly identified (FN=0, no missed examples).

[0202] The precision is:

[0203]

[0204] In the formula, Precision ranges from [0,1], where 0 indicates that all predicted positive examples are misclassified (TP=0, all predicted positive examples are negative examples), and 1 indicates that all predicted positive examples are correct (FP=0, no misclassified positive examples).

[0205] F1 value:

[0206]

[0207] In the formula, the F1 value represents the harmonic mean of the two indicators, Precision and Recall, which means that the two indicators are equally important and are more sensitive to imbalance between the two indicators than the arithmetic mean. The range of the F1 value is [0,1], where 0 means that at least one of Precision and Recall is 0, and 1 means that both Precision and Recall are 1 (TP=1, FP=0, FN=0).

[0208] All metrics range from 0 to 1, with values ​​closer to 1 indicating better model performance.

[0209] Accuracy focuses on overall correctness and is suitable for scenarios with balanced classes. When the data is imbalanced, such as when there are few samples of extreme rainstorm events, accuracy may become ineffective. For example, in a scenario with 95% negative examples, predicting all of them as negative examples can achieve 95% accuracy, but the recall rate is 0. In this case, recall, precision, and F1 score are more likely to reflect the model's ability to predict positive examples.

[0210] Recall and precision focus on the "capturing ability" and "predictive reliability" of positive examples, respectively, and are suitable for evaluating extreme events (few positive examples).

[0211] The F1 score forces a balance between recall and precision, is suitable for class imbalance scenarios, and can be used as a core evaluation metric for extreme rainstorm forecasting.

[0212] Given the high-risk nature of extreme rainstorms in hilly areas, the recall rate and F1 score should be given special attention to avoid missed cases while maintaining a comprehensive balance.

[0213] Example 2:

[0214] This embodiment is an application example of data from Embodiment 1:

[0215] Xingguo County, Ganzhou City, Jiangxi Province, is the study area. This area has a typical red soil hilly landform with an average slope of 25° to 35° and an average annual precipitation of 1500 to 1800 mm. Extreme rainstorms (daily precipitation ≥ 100 mm) occur 2 to 3 times a year on average. Historically, severe flash floods and debris flows have been caused by rainstorms.

[0216] Meteorological data collection: Hourly meteorological data (12 characteristics including precipitation intensity, temperature, humidity, and wind speed) from 2015 to 2024, from Xingguo County Meteorological Station and ERA5 reanalysis data from ECMWF.

[0217] Outlier handling: Detect outliers in precipitation intensity, such as a precipitation intensity of 120 mm / h at a certain time in July 2023, which exceeds the dynamic threshold, and replace it with the 3-hour sliding window average.

[0218] Topographic data collection and processing: Three features, namely slope, aspect and altitude, were extracted from the 30m resolution digital elevation model (DEM) and combined with land use type (forest, cultivated land and water) to construct a composite topographic index.

[0219] Label data: A total of 280 samples were labeled with extreme rainstorm events (daily precipitation ≥100mm), with positive examples accounting for 15%, indicating class imbalance.

[0220] Dataset Splitting and Simulation Data Augmentation: The dataset was split into a training dataset and a test dataset in a 7:3 ratio. The training dataset contained 196 samples, and the test dataset contained 84 samples. Then, 1000 sets of simulated extreme rainstorm events were generated based on historical data. The generalization ability of the model was expanded by randomly perturbing the precipitation intensity (±20%) and terrain parameters (slope ±5°).

[0221] Model building, training, and parameter optimization: Radial basis function (RBF) kernel is selected as the kernel function of SVM.

[0222] RBF kernel SVM parameter optimization: The grid search range is C = [0.1, 1, 10, 100], gamma = [0.001, 0.01, 0.1, 1]. The optimal parameters were obtained through training: C = 10, gamma = 0.01, and the F1 score reached 0.85 through 5-fold cross-validation.

[0223] Comparative experiments showed that the linear kernel SVM's performance was limited due to the significant nonlinear characteristics of the terrain, resulting in an F1 score of 0.78; the polynomial kernel SVM, due to the tendency of high-order polynomials to overfit, resulted in an F1 score of 0.80. The advantage of the RBF kernel SVM lies in its ability to improve the separability of complex terrain features through nonlinear mapping, thus demonstrating stronger generalization ability on simulated data.

[0224] The model training is completed according to the following criteria: the rate of change of the objective function value is <10% over 50 consecutive iterations. -6 The convergence condition is met; the cross-validation F1 value fluctuates by less than 1%, and the accuracy of the training dataset remains stable at 93%. Furthermore, considering practical business applications, in terms of real-time inference efficiency, single-sample prediction takes less than 0.1 seconds, meeting the real-time requirements of the early warning system.

[0225] Model evaluation results: In terms of performance metrics on the test dataset, accuracy = 0.91, recall = 0.88, precision = 0.86, F1 = 0.87; it significantly reduces the risk of missed cases and reduces the social cost of misclassification, and its overall performance is better than traditional models.

[0226] Confusion Matrix and Error Analysis: The confusion matrix is...

[0227]

[0228] Missed assessments (FN=6) often occur in deep mountainous areas with complex terrain, such as Meijiao Town in Xingguo County, and need to be optimized by combining satellite remote sensing data. Misjudgments (FP=5) are mainly due to short-duration heavy rainfall (<3 hours) not reaching extreme standards, which can be improved by extending the time window (e.g., 6 hours).

[0229] Compared with existing systems:

[0230] (1) In this embodiment, the evaluation index F1 = 0.87. Its advantage lies in its strong nonlinear feature extraction capability, which reduces the false negative rate by 40% compared with the traditional model. Moreover, the model size is <1MB, which can run in real time on edge devices.

[0231] (2) A certain water conservancy department model has an evaluation index F1 = 0.72. The model is based on linear regression, which has poor adaptability to complex terrain and a high misjudgment rate.

[0232] (3) A certain meteorological bureau’s Lingxi big model has an evaluation index F1 = 0.82, which is excellent in forecasting heavy rainfall during the flood season. However, it relies on cloud computing, requires high computing power, and is difficult and costly to deploy edge devices. It is not suitable for areas with weak networks.

[0233] Finally, it should be noted that the above is only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the present invention.

Claims

1. A method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling, characterized in that, Includes the following steps: Step S1, Data Collection and Preprocessing: Collect the ERA5 dataset from the European Centre for Medium-Range Weather Forecasts (ECMWF), measured data from meteorological stations within the study area, and DEM (Digital Elevation Model) data reflecting the topographic features of the study area; preprocess the collected data, including quality control and standardization. The quality control includes removing outliers and filling missing values ​​using a hybrid strategy of spatiotemporal coupling dynamic thresholding, multivariate joint testing, spatiotemporal density clustering verification, and physical mechanism inversion verification. The preprocessed data is then divided into training datasets and test datasets. The physical mechanism inversion verification involves measuring the water vapor flux of the outliers to be removed. Verify and determine whether it exceeds the 95th percentile of the historical extreme values ​​of the study area. If it does, it is determined to be a real extreme event, and the data is retained. Step S2, Multi-factor Feature Extraction: Principal component analysis is used for dimensionality reduction to extract key meteorological features and calculate the topographic index. To illustrate the influence of topography on precipitation: In the formula, a is the catchment area and b is the slope; Step S3, Key Influencing Factor Screening: The random forest algorithm is used to rank the importance of the extracted multi-factor features, and features with higher importance scores are selected as key influencing factors for extreme rainstorm events; including: Step S31, Calculate overall importance: Construct multiple decision trees using random forest to classify or regress the data. Each decision tree randomly selects a subset of samples and features for training. After training, the overall importance of each feature is evaluated by calculating the average reduction in impurity across all decision trees. Let feature... In the The importance of each decision tree is ,feature The formula for calculating the overall importance is: In the formula, Features The overall importance, For the number of decision trees, Features In the The importance of each decision tree; Step S32, Select key influencing factors: Sort the multi-factor features extracted in step S2 according to their overall importance from largest to smallest, and select the top ones with higher overall importance. Several characteristics serve as key influencing factors for extreme rainstorm events; among them... The range of values ​​is , The total number of multi-factor features extracted in step S2; the top features with higher overall importance are selected. When considering a feature, ensure that the cumulative importance is greater than or equal to the second preset threshold. The formula for cumulative importance is: In the formula, Features for sorting from largest to smallest The overall importance; Step S4, Prediction Model Construction and Training: A prediction model for extreme rainstorm events is constructed using the support vector machine algorithm; the key influencing factors of the selected extreme rainstorm events are used as the input features of the prediction model, and "whether the extreme rainstorm event has occurred" is used as the output label of the prediction model. The prediction model is trained until the model training is completed. Step S5, Evaluation of Prediction Results: The extreme rainstorm event prediction model trained in Step S4 is evaluated using the test dataset described in Step S1. The evaluation metrics include accuracy, recall, precision, and a comprehensive metric balancing recall and precision. value; The specific steps for removing outliers in step S1 include: The first step, univariate spatiotemporal dynamic threshold detection: construct a spatiotemporal sliding window for each meteorological variable and calculate the dynamic threshold. : In the formula, It is the 75th percentile within the spatiotemporal window. Interquartile range, , The 25th percentile of the spatiotemporal window in mainland China; This represents the average elevation difference between the station and its neighboring areas. For reference elevation parameters, As a terrain adjustment factor, For the station elevation, The average elevation of the study area, Standard deviation of elevation; For suspected outliers exceeding the threshold, the time series curvature is used for analysis. Verification, the formula is: A value is considered an outlier if it meets the following conditions: In the formula, The time step is the current moment. For meteorological variables in time Observations at that location For time The observation at the next time step, time The observation from the previous time step, For the standard deviation of curvature, This is the seasonal average. The second step, multivariate joint anomaly detection: calculating the multivariate Mahalanobis distance. : In the formula, For multivariate observation vectors, It is the mean vector. This is the transpose operator for matrices or vectors. It is the covariance matrix; Introducing terrain weight correction threshold : In the formula, For the number of variables, For degrees of freedom The 99th percentile of the chi-square distribution, This refers to the terrain index at the current location. To determine the maximum topographic index for the study area, For the catchment area, Slope; If satisfied This is then identified as a multivariate coupling outlier; The third step is to verify the spatiotemporal density clustering: calculate the composite spatiotemporal distance metric. : In the formula, For spatial Euclidean distance, For time intervals, To study the maximum Euclidean distance between all pairs of spatial points within the region, To study the maximum time span within a time range, The Pearson correlation coefficient is used to determine the trend of the variable. , , For adjustable weighting coefficients, ; Calculate the adaptive density parameters: In the formula, The neighborhood radius, The mean of the distance distribution. MinPts is the minimum number of cluster points, representing the standard deviation of the distance. This represents the total number of samples. If a data sample point does not belong to any density cluster and is not covered by more than 80% of its neighboring points, it is considered an outlier. Four steps, physical mechanism inversion verification: Water vapor flux is used to verify outliers to be removed. Verification, the formula is: If the following conditions are met, the event is considered a real extreme event, and the data will be retained: In the formula, The 95th percentile of the historical extreme values ​​for the study area. Specific humidity is the mass of water vapor in a unit mass of moist air. A three-dimensional wind speed vector. The direction is east-west. In the south-north direction, Vertical direction The wind speed is in the east-west direction, with east being the positive direction. The wind speed is in the north-south direction, with north being the positive direction. Vertical wind speed, upward is positive.

2. The method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling as described in claim 1, characterized in that, The missing value filling in step S1 specifically involves using linear interpolation to fill in the missing values. Let the original data sequence be... The linear interpolation formula is: In the formula, For missing values, for Previously adjacent valid data, for The next valid data, Location of missing values Up to the previous valid data The next adjacent valid data The number of intervals.

3. The method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling as described in claim 1, characterized in that, The standardization process described in step S1 specifically involves transforming the values ​​of each variable into a standard normal distribution with a mean of 0 and a standard deviation of 1, in order to eliminate the influence of different variable units. The standardization formula is as follows: In the formula, For the first The first sample The standardized values ​​of the variables For the first The first sample One variable value, For the first The mean of each variable, For the first The standard deviation of each variable.

4. The method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling as described in claim 1, characterized in that, The multi-factor feature extraction in step S2 specifically includes the following sub-steps: Step S21, calculate the covariance matrix S: Let the measured meteorological station data of the training dataset mentioned in step S1 be... Data matrix ,in For the sample size, The total number of meteorological features; calculated Covariance matrix S: In the formula, This is the transpose operator for matrices or vectors. for The transpose matrix and covariance matrix for A real symmetric matrix; Step S22, Eigenvalue decomposition: The covariance matrix obtained in step S21... There exist orthogonal matrices and diagonal matrix , so that: Among them, the diagonal matrix The diagonal elements are the covariance matrix. eigenvalues orthogonal matrix The column vectors are the corresponding unit orthogonal eigenvectors. ; Step S23, Dimensionality Reduction: Select the previous... The eigenvectors corresponding to each eigenvalue They are used as column vectors to form the projection matrix. : in, The range of values ​​is Before selection When there are 10 eigenvalues, ensure that the cumulative variance contribution ratio is greater than or equal to the first preset threshold. The formula for the cumulative variance contribution ratio is: In the formula, Let i be the i-th eigenvalue sorted from largest to smallest; Through projection matrix The matrix of meteorological station measured data from the training dataset described in step S1 Projecting onto a lower-dimensional space yields the principal component matrix Y: Finally, the top k main meteorological features were extracted; Step S24, Calculate the terrain index: Calculate the terrain index using the digital elevation model data of the training dataset described in step S1. .

5. The method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling as described in claim 1, characterized in that, The support vector machine algorithm described in step S4 is as follows: For the linearly separable case, the formula for solving the optimization problem using support vector machine is: In the formula, w is the weight vector, and b is the bias term. For the input sample, A function that maps input samples to a high-dimensional feature space. For the category label of the sample, ; For the linearly inseparable case, the formula for solving the optimization problem using support vector machines is: In the formula, As slack variables, This is a penalty factor.

6. The method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling as described in claim 1, characterized in that, The specific process of training the prediction model until the model training is completed in step S4 is as follows: During the training process, grid search is used to traverse the kernel function parameter combinations, and cross-validation is used to optimize the model parameters and select the optimal penalty factor. and kernel function parameters; The completion of model training specifically includes the following criteria: The optimization algorithm meets the convergence conditions: the change in the objective function between adjacent iterations is less than a third preset threshold; the training samples satisfy the Karush-Kuhn-Tucker conditions. Preset upper limit for the number of iterations: If the optimization algorithm fails to converge within a fixed number of iterations, it will be forcibly stopped to avoid infinite loops; Stable performance on the test dataset: If the model performance no longer improves after a certain number of consecutive rounds on the cross-validation or the test dataset, training will be terminated early to prevent overfitting. Parameter optimization complete: The optimal penalty factor was found through grid search. Training is considered complete when the model's solution on the training set satisfies the constraints of the optimization problem, along with the kernel function parameters.

7. The method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling as described in claim 1, characterized in that, The formula for calculating accuracy in step S5 is as follows: In the formula, For the number of true cases, The number of true counterexamples The number of false positives. The number of false counterexamples; The range is 0 indicates that the prediction of the extreme rainstorm event was completely wrong, and 1 indicates that the prediction of the extreme rainstorm event was completely correct.

8. The method for predicting and analyzing extreme rainstorm events in hilly areas based on multi-factor coupling as described in claim 1, characterized in that, The F1 value mentioned in step S5 is: In the formula, Value representation and The harmonic mean of the two indicators means that both indicators are equally important; The range of values ​​is 0 represents and At least one is 0, 1 indicates and All are 1.

Images