A rainfall prediction method based on signal stationary and non-stationary feature separation double model fusion

Abstract

The application discloses a rainfall prediction method based on signal stationary and non-stationary feature separation double-model fusion, first acquires and pre-processes rainfall prediction related data; then separates stationary feature sequence and non-stationary feature sequence from the pre-processed Beidou original observation data through a singular spectrum analysis method; then constructs a first prediction model and a second prediction model, weights and fuses the output stationary feature prediction result and the non-stationary feature prediction result to obtain the final prediction value of a traditional prediction variable; then acquires a rainfall prediction factor and a rainfall type determination factor; finally constructs and trains a comprehensive rainfall prediction model, takes the final prediction value of the traditional prediction variable, the receiver hardware differential code bias time sequence and the water vapor vertical distribution feature as a prediction feature set, and then outputs a rainfall prediction result through the pre-trained comprehensive rainfall prediction model.

Description

Technical Field

[0001] This invention relates to the interdisciplinary field of navigation satellite remote sensing and weather forecasting, specifically a dual-model fusion forecasting method based on the separation of stationary and non-stationary features of signals. Background Technology

[0002] Rainfall plays a crucial role in atmospheric water cycle and meteorological research. Accurate and timely rainfall forecasts can mitigate the damage caused by floods and droughts and ensure public safety during severe weather. GNSS (Global Navigation Satellite System) signals undergo refraction and delay during propagation in the troposphere; this delay is known as tropospheric delay. Tropospheric delay primarily refers to ZTD (Zenith Total Delay), which can be further divided into ZHD (Zenith Hydrostatic Delay) and ZWD (Zenith Wet Delay). Zenith hydrostatic delay is caused by dry gases such as oxygen and nitrogen, while zenith wet delay is related to water vapor content and can be used to derive PWV (Precipitable Water Vapor). Ground-based GNSS, with its low cost, all-weather availability, high temporal resolution, and high accuracy, has become an important means of obtaining atmospheric water vapor information. Utilizing atmospheric parameters such as ZTD, ZWD, and PWV obtained through GNSS inversion for rainfall forecasting is currently a research hotspot in the field of GNSS meteorology.

[0003] Existing GNSS rainfall prediction methods can be broadly categorized into four types. The first type is based on thresholds or empirical values. It uses statistical analysis of the correlation between atmospheric parameters retrieved from GNSS and rainfall events to set thresholds or construct empirical indices for rainfall prediction. This type of method is computationally simple and has clear physical meaning, but it does not consider the non-stationary characteristics of the signal, resulting in low prediction accuracy. The second type is based on machine learning or deep learning. It employs machine learning models such as support vector machines, random forests, and neural networks, using time series data of GNSS meteorological parameters as input features for rainfall prediction. This type of method has strong nonlinear fitting capabilities, but it only eliminates receiver hardware bias (DCB) as an error, without incorporating DCB as a predictive factor into the model.

[0004] The third type of method is based on Numerical Weather Prediction (NWP) data assimilation. This involves assimilating tropospheric data retrieved from GNSS and incorporating it into the numerical weather prediction model, thereby improving the accuracy of rainfall forecasts. This type of method has strong physical consistency and continuous spatial coverage, but the assimilated water vapor information is mainly based on whole-layer integration, making it difficult to accurately characterize changes in local small-scale water vapor structure. The fourth type of method is based on multi-source data fusion.

[0005] This method integrates GNSS data with various sources, including satellite remote sensing, weather radar, and ground meteorological observations, to construct a more comprehensive rainfall prediction model. While this approach leverages the complementary advantages of multi-source data, it only considers the spatial and temporal fluctuations of the predictor variables, neglecting the impact of altitude distribution on rainfall. Among the four types of GNSS rainfall prediction algorithms mentioned above, only receiver hardware bias (DCB) is treated as an error to be eliminated. However, research shows that when air humidity changes significantly, the receiver hardware bias DCB fluctuates, and this fluctuation can be considered a precursor to rainfall. Furthermore, the time series of atmospheric parameters retrieved from GNSS exhibits non-stationary characteristics, but none of these four rainfall prediction algorithms consider these non-stationary features or perform targeted separation processing, resulting in low rainfall prediction accuracy. In addition, these four types of GNSS rainfall prediction algorithms only consider the spatial and temporal fluctuations of the predictor variables, neglecting the impact of altitude distribution on rainfall. However, existing research shows that the altitude of water vapor distribution directly determines the type of future rainfall. Water vapor concentrated in the mid-to-high atmosphere easily forms strong convection, leading to extreme rainfall, while water vapor concentrated near the ground is mostly stable precipitation. Therefore, researching a BeiDou satellite rainfall prediction algorithm that can separate stationary and non-stationary characteristics of the predicted quantity, use DCB bias as a prediction factor, and introduce the vertical characteristic distribution of water vapor can effectively improve the accuracy and timeliness of rainfall prediction, and has important theoretical significance and value. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of existing technologies and propose a rainfall prediction method based on dual-model fusion of stationary and non-stationary signal features. This method can effectively overcome the problems of insufficient accuracy and inadequate information utilization in GNSS rainfall prediction methods. The fused receiver hardware bias DCB is based on the rainfall precursor information, using the DCB fluctuation as a reference. The DCB bias, along with traditional predictor variables such as ZTD, ZWD, and PWV, are used as input features to construct a rainfall prediction model, improving prediction accuracy. The dual-model fusion prediction algorithm based on the separation of stationary and non-stationary signal features first uses a time-series decomposition method to separate stationary and non-stationary features from the time series of traditional predictor variables such as ZTD, ZWD, and PWV. The separated stationary feature signals are then predicted using a mixed linear prediction model (MixLinear is an extremely low-resource multivariate time series prediction model proposed at ICLR 2026). The separated non-stationary signal features are then predicted using a destationary Fourier transform and coefficient network prediction model (DFCNformer).

[0007] The De-stationary Fourier and Coefficient Network Transformer (a Transformer framework for non-stationary time series forecasting) prediction algorithm improves prediction accuracy. Introducing the vertical distribution characteristics of water vapor directly and effectively determines the type of future rainfall, further enhancing the prediction accuracy of extreme rainfall. The fusion of these three methods significantly improves the reliability, accuracy, and timeliness of BeiDou satellite rainfall prediction.

[0008] To achieve the above objectives, the technical solution specifically adopted by the present invention is as follows:

[0009] A rainfall prediction method based on the fusion of dual models that separate stationary and non-stationary signal features includes the following steps:

[0010] Step 1: Acquire and preprocess rainfall forecasting-related data, which includes BeiDou raw observation data, historical radiosonde data, and reanalysis datasets;

[0011] Step 2: Separate stationary and non-stationary feature sequences from the preprocessed BeiDou raw observation data using singular spectrum analysis.

[0012] Step 3: Construct a first prediction model to predict the stationary feature sequence and obtain the stationary feature prediction result; construct a second prediction model to predict the non-stationary feature sequence and obtain the non-stationary feature prediction result; weight and fuse the stationary feature prediction result and the non-stationary feature prediction result to obtain the final predicted value of the traditional predictor variable.

[0013] Step 4: Extract the receiver hardware differential code deviation time series from the BeiDou observation data, and use it as a rainfall prediction factor after preprocessing.

[0014] Step 5: Obtain the vertical distribution data of atmospheric precipitable water from the reanalysis dataset, calculate the vertical distribution data of zenith wet delay and zenith tropospheric delay from the radiosonde data, and extract the vertical distribution characteristics of water vapor as a precipitation type determination factor.

[0015] Step 6: Construct and train a comprehensive rainfall prediction model. Use the final predicted values ​​of traditional prediction variables, the receiver hardware differential code deviation time series, and the water vapor vertical distribution characteristics as the prediction feature set, and then output the rainfall prediction result through the pre-trained comprehensive rainfall prediction model.

[0016] Preferably, the specific steps for acquiring and preprocessing rainfall forecast data include:

[0017] The original BeiDou observation data was processed using precise point positioning technology, the cutoff elevation angle was set, the ionospheric delay effect was eliminated by using an ionospheric-free combination model, and the zenith tropospheric delay was estimated by using broadcast ephemeris files and precise orbit and clock difference files through Kalman filtering.

[0018] Based on the air pressure, latitude, and altitude at the receiver site, the zenith static delay is calculated using the Sastamonin model. The zenith wet delay is obtained by subtracting the zenith static delay from the zenith tropospheric delay. Precipitation is calculated based on the zenith wet delay and the conversion factor. The zenith tropospheric delay, zenith wet delay, and precipitation constitute traditional rainfall prediction variables.

[0019] Vertical profile data of temperature, humidity and air pressure are extracted from the sounding database, and atmospheric precipitable water and pressure layer data are obtained from the reanalysis dataset to obtain water vapor vertical distribution data.

[0020] Preferably, in step 2, the preprocessed zenith tropospheric delay time series, zenith wet delay time series, and precipitable water time series are decomposed into stationary and non-stationary characteristic series.

[0021] Preferably, step 2, the specific steps for separating stationary and non-stationary feature sequences, include:

[0022] Singular spectral analysis is used for feature separation. The window length is set, and the time series is mapped to a trajectory matrix by an embedding operator.

[0023] Singular value decomposition is performed on the trajectory matrix to obtain multiple components;

[0024] Components associated with steady trends and periodicity are grouped into one group, while components associated with noise and non-stationary abrupt changes are grouped into another group.

[0025] The two components are reconstructed into one-dimensional time series using the diagonal averaging method, yielding stationary and non-stationary characteristic sequences.

[0026] Preferably, the specific steps for constructing a first prediction model to predict stationary feature sequences and obtaining stationary feature prediction results include:

[0027] The first prediction model is a hybrid linear prediction model, which includes a time-domain segmented processing path and a frequency-domain adaptive low-rank filtering path.

[0028] The stationary feature sequence is input into the hybrid linear prediction model. Local waveform features are extracted and cross-segment dependencies are captured through the time-domain segmentation processing path. The downsampled data sequence is then subjected to Fourier transform and low-rank decomposition adaptive spectral filtering through the frequency-domain adaptive low-rank filtering path.

[0029] The outputs of the two paths are upsampled and then fused by addition to obtain the prediction results of the stationary feature components.

[0030] As a preferred embodiment, the extraction of local waveform features and capture of cross-segment dependencies through time-domain segmentation processing specifically includes:

[0031] The input sequence is downsampled and divided into multiple non-overlapping time segments;

[0032] Perform intra-segment linear projection on each time segment to extract local waveform features;

[0033] After concatenating all intra-segment features, perform inter-segment linear projection to capture cross-segment dependencies.

[0034] Preferably, the Fourier transform and low-rank decomposition adaptive spectral filtering of the downsampled data sequence via the frequency domain adaptive low-rank filtering path specifically includes:

[0035] Apply the Fast Fourier Transform to the downsampled data sequence to obtain the Fourier Transform result;

[0036] Adaptive spectral filtering is achieved by employing low-rank decomposition technology. The spectrum of each frequency band is mapped to a low-dimensional latent space and reconstructed to obtain the frequency domain representation after low-rank adaptive spectral filtering.

[0037] Preferably, the specific steps for constructing a second prediction model to predict the non-stationary feature sequence and obtaining the non-stationary feature prediction result include:

[0038] The second prediction model is a destationary Fourier transform and coefficient network prediction model;

[0039] The non-stationary feature sequences are normalized and stabilized to ensure that the statistical properties of the input data remain consistent.

[0040] The processed sequences were decomposed into seasonal and trend components;

[0041] For the seasonal components, a de-stationary Fourier attention mechanism is used, and an encoder-decoder architecture is employed for prediction to obtain the prediction results of the seasonal components.

[0042] The trend components are predicted using a multilayer perceptron and a dual-coefficient network to obtain the prediction results.

[0043] The prediction results of seasonal components are added together with the prediction results of trend components, and then inverse stationarization is performed to obtain the prediction results of non-stationary characteristic components.

[0044] Preferably, the processed sequence is decomposed into seasonal and trend components, specifically including:

[0045] By using moving average filters with different window sizes to capture various trend patterns, and by using adaptive weights to merge all trend patterns, the final trend component is obtained.

[0046] The final seasonal component is obtained by subtracting the final trend component from the original processed sequence.

[0047] Preferably, the specific steps for predicting the seasonal components using a de-stationary Fourier attention mechanism and an encoder-decoder architecture include:

[0048] The seasonal components are processed through a linear layer to obtain a query matrix, a key matrix, and a value matrix.

[0049] Calculate the non-stationary factor, which includes the standard deviation of the input sequence and the mean value of the query matrix over the time dimension;

[0050] The query matrix, key matrix, and value matrix are transformed into the frequency domain using Fourier transform;

[0051] An attention mechanism is applied in the frequency domain, and the non-stationary factor is used to correct the attention calculation results;

[0052] The frequency domain results are converted back to the time domain by inverse Fourier transform, and then passed to the multilayer encoder and multilayer decoder for processing to obtain the prediction results of the seasonal components.

[0053] Preferably, the specific steps for predicting the trend component using a multilayer perceptron and a dual-coefficient network include:

[0054] The trend component sequence is passed through the input coefficient network to obtain the horizontal parameters and fluctuation parameters of the input data, and then the output is normalized.

[0055] The normalized output is input into the multilayer perceptron to complete model training;

[0056] After the multilayer perceptron outputs the prediction result, it predicts future values ​​through the output parameter network and outputs the normalized prediction result.

[0057] The normalized prediction results are then subjected to inverse normalization to obtain the prediction results of the trend components.

[0058] Preferably, the specific steps for weighted fusion of the stationary feature prediction results and the non-stationary feature prediction results to obtain the final predicted value of the traditional predictor variable include:

[0059] The weights of stationary and non-stationary feature components are determined by using Bayesian optimization methods to minimize the prediction error of the validation set.

[0060] Based on the weights of the stationary feature components and the non-stationary feature components, the prediction results of the stationary feature and the non-stationary feature are weighted and summed to obtain the final predicted value of the traditional predictor variable, wherein the sum of the weights of the stationary feature components and the non-stationary feature components is 1.

[0061] As a preferred method, the specific steps for extracting the receiver hardware differential code offset time series from BeiDou observation data and preprocessing it as a rainfall prediction factor include:

[0062] Based on the dual-frequency ranging code and carrier phase observations, the differential code deviation information is extracted using the carrier phase smoothing pseudorange method to eliminate ionospheric delay and obtain the observed value of the sum of satellite differential code deviation and receiver differential code deviation.

[0063] Add a zero-mean centroid reference and separate the differential code deviation time series of each receiver by least squares adjustment;

[0064] Outliers were removed using statistical methods, missing values ​​were imputed using linear interpolation, and high-frequency noise was removed using moving average filtering to obtain the preprocessed receiver hardware differential code deviation time series.

[0065] The receiver hardware differential code bias time series is decomposed into stationary and non-stationary feature components, which are then concatenated with the stationary and non-stationary feature sequences of the traditional rainfall prediction variables, respectively, and used as input to the prediction model.

[0066] Preferably, the specific steps for obtaining the vertical distribution data of atmospheric precipitable water from the reanalysis dataset, calculating the vertical distribution data of zenith wet delay and zenith tropospheric delay from radiosonde data, and extracting water vapor vertical distribution characteristics as a precipitation type determination factor include:

[0067] Download the specific humidity, relative humidity, temperature and gravitational potential variables of different pressure layers from the reanalysis dataset to obtain the vertical distribution data of atmospheric precipitable water.

[0068] Temperature, air pressure and specific humidity of each altitude layer are extracted from the sounding database, and the zenith wet delay and zenith tropospheric delay of each altitude layer are calculated by integration.

[0069] After unifying the elevation datum, a tomographic grid is defined. In the horizontal direction, bilinear interpolation is used to interpolate the data to the horizontal grid points of the target grid, and in the vertical direction, exponential interpolation is used to resample the data to the target standard height layer.

[0070] The vertical distribution characteristics of water vapor are calculated based on the interpolated data. These characteristics include the vertical gradient of water vapor, the height center of the vertical distribution of water vapor, the proportion of water vapor in the lower layer, the proportion of water vapor in the middle layer, and the proportion of water vapor in the upper layer.

[0071] As a preferred option, the specific steps for constructing a comprehensive rainfall prediction model include:

[0072] The samples are divided into training set, validation set and test set in chronological order, and the historical time window length, sampling interval and prediction step size are set.

[0073] Construct a comprehensive rainfall prediction model that combines classification and regression models;

[0074] For the classification model, the vertical distribution characteristics of water vapor are used as input variables. Historical extreme rainfall events are labeled as specified categories according to the set extreme rainfall threshold. A random forest algorithm is used to construct a rainfall discrimination model, and the final predicted category is determined by majority voting.

[0075] For the regression model, the traditional rainfall prediction variables and the receiver hardware differential code deviation time series are used as input variables. The first prediction model is used to predict stationary series, and the second prediction model is used to predict non-stationary series. The prediction results are then weighted and fused to obtain the final rainfall prediction result.

[0076] The discrimination result of the classification model and the prediction result of the regression model are used together as the output of the integrated rainfall prediction model.

[0077] Preferably, the random forest algorithm generates a training subset by sampling with replacement and randomly selects a feature subset at each node for splitting, thereby reducing the correlation between trees and improving generalization performance.

[0078] This invention has the following characteristics and beneficial effects:

[0079] This invention provides a BeiDou satellite precipitation prediction method that integrates receiver hardware bias DCB, a dual-model fusion prediction algorithm based on the separation of stationary and non-stationary signal features, and incorporates water vapor vertical distribution characteristics. Compared with existing technologies, it has the following advantages:

[0080] First, this invention effectively solves the problem of low rainfall prediction accuracy caused by neglecting the non-stationary features of signals in existing methods by separating stationary and non-stationary features of traditional predictive variables such as ZTD, ZWD, and PWV, and then using the MixLinear model and DFCNformer model for targeted prediction. The MixLinear model can efficiently capture local time-series features and global trends in stationary signals, while the DFCNformer model effectively addresses the inherent non-stationarity of non-stationary time series through a destationary Fourier attention mechanism and a dual-coefficient network. The weighted fusion of the two models can comprehensively improve the accuracy and stability of rainfall prediction.

[0081] Second, this invention incorporates receiver hardware bias (DCB), traditionally considered an error source and eliminated, into the rainfall prediction model as a rainfall prediction factor. By extracting the time-series features of DCB and separating them into stationary and non-stationary types, and then concatenating them with their corresponding feature components before inputting them into the prediction model, it is equivalent to adding an independent information source. This fully utilizes the fluctuations in DCB caused by significant changes in air humidity—a precursor to rainfall—effectively solving the problems of insufficient information utilization and lack of precursory rainfall information in existing methods, and significantly improving the model's ability to detect rainfall events in advance.

[0082] Third, this invention systematically introduces the vertical distribution characteristics of water vapor as a factor for determining rainfall type. By obtaining the vertical distribution data of PWV from ERA5 reanalysis data and calculating the vertical distribution data of ZWD and ZTD from radiosonde data, features such as the vertical gradient of water vapor, the center height of the vertical distribution of water vapor, the proportion of water vapor in the lower, middle, and upper layers are extracted. This enables the prediction model to better distinguish between extreme rainfall and stable precipitation formed by strong convection, effectively solving the problem that existing methods ignore the influence of water vapor height distribution on rainfall, significantly improving the ability to identify and perceive extreme rainfall events, and reducing the false alarm rate of extreme rainfall.

[0083] Fourth, this invention adopts a comprehensive rainfall prediction framework that combines classification and regression models. The classification model uses the vertical distribution characteristics of water vapor to determine whether it is extreme rainfall, while the regression model uses a fusion model of stationary and non-stationary rainfall to predict rainfall. The two work together to accurately distinguish rainfall types and output quantitative rainfall prediction results, providing more efficient and robust technical support for short-term weather forecasting, extreme weather warnings, disaster prevention and mitigation, and agricultural production. It has broad application prospects and significant socio-economic value. Attached Figure Description

[0084] Figure 1 This is a flowchart of the data acquisition and preprocessing process for rainfall prediction in this embodiment.

[0085] Figure 2 This is a flowchart of a dual-model fusion prediction method based on the separation of stationary and non-stationary features of a signal in this embodiment.

[0086] Figure 3 This is a flowchart of the rainfall prediction variable prediction process based on DCB feature fusion in this embodiment.

[0087] Figure 4 This is a flowchart of the extreme rainfall discrimination and prediction process based on the vertical distribution characteristics of water vapor in this embodiment. Detailed Implementation

[0088] The present invention will now be described in detail with reference to specific embodiments. These embodiments will help those skilled in the art to further understand the present invention, but do not limit the invention in any way. It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0089] A rainfall prediction method based on dual-model fusion of stationary and non-stationary signal features is proposed. This method integrates receiver hardware bias (DCB), a dual-model fusion prediction algorithm based on stationary and non-stationary signal features, and a BeiDou satellite rainfall prediction method incorporating water vapor vertical distribution characteristics. This effectively improves the reliability, accuracy, and timeliness of BeiDou satellite rainfall prediction. By separating stationary and non-stationary features from traditional rainfall prediction variables, the MixLinear prediction algorithm is used for stationary features, while the DFCNformer prediction algorithm is used for non-stationary features. This effectively addresses the problem of low rainfall prediction accuracy caused by traditional algorithms not considering non-stationary features. The MixLinear algorithm organically combines segmented trend extraction in the time domain with adaptive low-rank spectral filtering in the frequency domain, effectively capturing local time-series features and achieving high computational efficiency while maintaining prediction accuracy. The DFCNformer model uses a non-stationary Fourier attention mechanism for seasonal time series and a dual-coefficient network and multilayer perceptron to predict long-term trends for trend time series. Combining these two approaches with stationarization effectively solves the inherent non-stationarity problem of non-stationary time series. By fusing receiver hardware bias DCB (Distributed Direct Bias), the DCB, traditionally considered an error to be eliminated, is transformed into a precursor to rainfall, solving the problems of insufficient information utilization and lack of rainfall precursor information. By introducing water vapor vertical distribution characteristics, using ERA5 reanalysis data to obtain PWV (Pulse Wave Volume) vertical distribution data, and using radio sounding (RS) data interpolation to obtain ZWD (Zero Distance) and ZTD (Zero Distance Tolerance) vertical distribution data, the prediction model can better perceive the distribution of water vapor at different altitudes, thereby better distinguishing between extreme rainfall and stable precipitation formed by strong convection, effectively differentiating rainfall types, and improving the accuracy of extreme rainfall prediction. This invention effectively improves the accuracy of existing GNSS rainfall prediction by separating and predicting traditional rainfall prediction variables, incorporating the time-varying characteristics of receiver hardware bias DCB and water vapor vertical distribution characteristics into the prediction factors, fully exploring rainfall precursor information, accurately distinguishing rainfall types, and providing reliable technical support for short-term rainfall forecasting, extreme weather warnings, and disaster prevention and mitigation. The steps include:

[0090] Step 1: Acquire and preprocess rainfall prediction-related data, which includes BeiDou raw observation data, historical radiosonde data, and reanalysis datasets.

[0091] In this embodiment, as Figure 1 As shown, precise point positioning (PPP) calculations are first performed using the open-source processing software RTKLIB, i.e., a 5° cutoff elevation angle is set, and an ionospheric-free combined model is used to eliminate the influence of ionospheric delay. Then, broadcast ephemeris files and precise orbit and clock error files provided by IGS are applied. Finally, Kalman filtering is used to estimate... . It can be further decomposed into and The specific formula is as follows:

[0092]

[0093] In the formula, For zenith tropospheric delay ( ), For zenith wet delay ( ), Zenith static delay ( ), and so on , and The meanings of all are the same as here. Among them, The Saastamoinen model can be used for calculation; the specific formula is as follows:

[0094]

[0095] In the formula, Atmospheric pressure at the Beidou receiver ( ), Latitude of the receiver site For the receiver altitude ( ).but can be minus , It is possible The specific formula is as follows:

[0096]

[0097]

[0098] In the formula, It is a conversion factor (dimensionless, It is a transformation factor (dimensionless)), and formulas (3) and (4) are... They mean the same thing. The density constant of liquid water ( ),Pick , The specific gas constant of water vapor ( ),Pick , and Let be an atmospheric physical constant, and take . ,Pick , Water vapor weighted average temperature ( ). , , These constitute the traditional variables for rainfall prediction. Data for regional stations were downloaded from the University of Wyoming Sounding Database website, and vertical profile data such as temperature, humidity, and air pressure were extracted. Data was also downloaded from ECMWF fifth-generation global climate reanalysis data. By extracting, analyzing, and performing horizontal and vertical compensation on barometric data, water vapor vertical distribution data can be obtained. Based on historical observation data, the historical time window length is set to 6 hours, the sampling interval to 5 minutes, and the prediction step size is set to 1 hour, 2 hours, and 3 hours respectively. The samples are then divided into a 60% training set, a 20% validation set, and a 20% test set in chronological order. Then... , , Time series data are Z-Score standardized (standard score standardization).

[0099] Step 2: Separate stationary and non-stationary feature sequences from the preprocessed BeiDou raw observation data using singular spectrum analysis.

[0100] In this embodiment, the predictor variable is decomposed into stationary and non-stationary feature sequences using a time-series decomposition method. Singular Spectrum Analysis (SSA) is employed for feature separation; SSA is a non-parametric and adaptive time-series decomposition method. The SSA process involves first setting... Given a time series of length N, For the window length, take The trajectory matrix is ​​constructed from the embedding operator, which can map the time series to... OK The formula for the Henkel matrix is ​​as follows:

[0101]

[0102]

[0103] In the formula, The trajectory matrix constructed for the embedding operator, ,..., N represents the original sequence observations, L is the window length for SSA decomposition, K is the number of columns in the trajectory matrix, and N is the total length of the time series. The meanings of N, L, and K in formulas (5) and (6) are the same. Next, Singular Value Decomposition (SVD) is performed, with the following formula:

[0104]

[0105] In the formula, The trajectory matrix is ​​the same as that in formula (6), where d is the number of non-zero singular values. For matrix The orthogonal system formed by the left singular vectors, For matrix An orthogonal system consisting of right singular vectors. The summation is the square of the singular values, and the superscript T signifies matrix transpose. Then, the components obtained from the SVD decomposition are divided into two groups: one group consists of components with the same stationary trend and those correlated with periodicity, and the other group consists of components correlated with noise and non-stationary abrupt changes. Finally, both groups are reconstructed into one-dimensional time series using the diagonal averaging method to obtain the stationary feature series. Non-stationary feature sequences .

[0106] Step 3: Construct a first prediction model to predict the stationary feature sequence and obtain the stationary feature prediction result; construct a second prediction model to predict the non-stationary feature sequence and obtain the non-stationary feature prediction result; weight and fuse the stationary feature prediction result and the non-stationary feature prediction result to obtain the final predicted value of the traditional predictor variable.

[0107] The first prediction model is a mixed linear prediction model (MixLinear model); the second prediction model is a destationary Fourier transform and coefficient network prediction model (DFCNformer model).

[0108] In this embodiment, as Figure 2 As shown, the MixLinear model is a bi-AND architecture that handles local trends in the time domain through linear transformation decomposition and global trends in the frequency domain through adaptive low-rank spectral filtering. To construct the MixLinear model, the input sequence is first... , This is the length of the historical window, and subsequent... The meaning is the same as here. Given the feature dimension, the MixLinear model generates prediction results via dual paths in the time and frequency domains, as shown in the following formula:

[0109]

[0110] In the formula, For time-domain segmentation processing path, This is a frequency-domain adaptive low-rank filtering path. Given a stationary feature sequence as input, For output, , To predict the step size, For feature dimensions, subsequent The meanings are the same as those here.

[0111] In the time domain, the input is downsampled by coefficients. Afterwards, it was divided into A length of The non-overlapping segments are defined by the following formula:

[0112]

[0113] In the formula, This is a collection of segments after segmentation. For the first A time segment, The length of the segment. This is the segment length, followed by... and The meanings are the same as here. Then, perform intra-segment linear projection on each segment to extract local waveform features. After concatenating the intra-segment features, perform inter-segment linear projection to capture cross-segment dependencies. The specific formula is as follows:

[0114]

[0115]

[0116]

[0117] In the formula, For the first Intra-segment linear projection of a time segment It is a linear projection function within the segment. Let be the feature dimension of the intrasegment projection. The result is the concatenation of features within all segments. This is the result of linear projection between segments. This is the inter-segment linear projection function.

[0118] In the frequency domain, a Fast Fourier Transform is applied to the downsampled data sequence, with the specific formula as follows:

[0119]

[0120] In the formula, The result is the Fourier transform. For the Fast Fourier Transform operator, For input The sequence was downsampled. Then, low-rank decomposition was used to implement adaptive spectral filtering, mapping each frequency band spectrum to the same low-dimensional latent space and reconstructing it. The specific formula is as follows:

[0121]

[0122] In the formula, The left matrix of the low-rank decomposition. The right matrix of low-rank decomposition, , Let be the dimension of the low-rank potential space. For complex fields, With formula ( In ) Same meaning This is the frequency domain representation after low-rank adaptive spectral filtering.

[0123] Finally, the outputs of the two paths are upsampled to restore them to the predicted length, and then fused by addition to obtain the final predicted value. The prediction result of the stationary feature component is output. .

[0124] Furthermore, to construct the DFCNformer model, firstly, the non-stationary feature sequences are... Normalization and stabilization are performed to ensure that the statistical characteristics of the input data remain consistent. The specific formula is as follows:

[0125]

[0126]

[0127]

[0128] In the formula, The mean vector of the input non-stationary feature sequence. Given the length of the input sequence, the subsequent... The meanings are the same as here. For the first A vector of observations at each time point For the standard deviation vector, the subsequent... The meanings are the same as here. The first normalized and stabilized version The values ​​at each time point are then decomposed into independent seasonal and trend components. Moving average filters with different window sizes are used to capture various trend patterns. Adaptive weights are then used to fuse all trend patterns to obtain the final trend component. The final seasonal component is obtained by subtracting the final trend component from the original time series. The specific formula is as follows:

[0129]

[0130]

[0131] In the formula, It is a seasonal ingredient. As a trend component, The input non-stationary feature sequence , For softmax operation, An adaptive weighting function related to the data. This is a moving average filter.

[0132] For seasonal components, a de-stationary Fourier Attention (DSF) mechanism is employed, utilizing a Transformer encoder-decoder architecture. The seasonal components are transformed into the frequency domain using a Fourier transform, and the attention mechanism is applied in the frequency domain to capture the importance of different frequency components. The specific formula for the attention mechanism is as follows:

[0133]

[0134]

[0135]

[0136]

[0137] In the formula, For standard attention mechanism functions, For the attention head dimension, the formula ( In ) With formula ( The meaning is consistent with that of ) These represent the query matrix, key matrix, and value matrix, respectively. This is the matrix transpose operator. and It is a non-stationary factor. In the time dimension The average value, seasonal components Obtained through a linear layer , For Fourier transform operators, This is the inverse Fourier transform operator. The output is then passed sequentially to the M-layer encoder and N-layer decoder for processing to obtain the prediction results of the seasonal components. .

[0138] For trend components, a Multi-Layer Perceptron (MLP) and a Dual Coefficient Network (Dual-CONET) are used for prediction. First, the stationary time series is passed through the BackConet input coefficient network, and the input sequence... and The learnable weights of the two coefficient networks within the l-th fully connected layer are defined by the following formula:

[0139]

[0140] In the formula, For the input coefficient network function, For the first The input sequence of a sample, Input window length, For the first A sample from time [time] At the time A continuous observation vector, and These represent the lateral parameters and fluctuation parameters of the input data in the input parameter network, respectively. For the first Learnable weights of fully connected layers For dimension function, For expectation operator, For the first Each sample at time... The observed values, For normalized output, and For learnable parameters, The total number of time series samples. For time step.

[0141] Then, the obtained Normalization is applied before inputting the data into the MLP to complete model training, thereby improving the MLP's prediction performance. Finally, the MLP outputs the prediction results. Then, the output parameter network is used to predict future values ​​and outputs the normalized prediction results. The specific formula for the output parameter network (PredConet) is as follows:

[0142]

[0143] In the formula, For the output parameter network function, and These are the level coefficients and amplitude coefficients predicted by the output parameter network, respectively. , , , , , , and The meanings are the same as those of formula (24). These are learnable weight parameters. For the first A sample from time [time] At the time Extended input sequence, For normalized output, The result was obtained through inverse normalization. Dual-CONET modules were added before and after the MLP. The overall trend prediction formula is as follows:

[0144]

[0145] In the formula, With formula (19) The meaning is the same. It is a two-coefficient network function. For multilayer perceptron functions, The input feature sequence is the trend component.

[0146] The seasonal and trend forecasts are added together, and then the results are processed through inverse stationarity to obtain the forecast results for the non-stationary characteristic components. At the same time, the DFCNformer model was constructed.

[0147] Finally, the stationary and non-stationary prediction results are weighted and fused to obtain the final prediction result of the traditional prediction quantity. Bayesian optimization is used to determine the weights, minimizing the prediction error on the validation set, and then the prediction results of the stationary feature components are combined. Prediction results of non-stationary eigencomponents The final predicted value is obtained by weighted fusion, and the specific formula is as follows:

[0148]

[0149]

[0150] In the formula, The weights of the stationary eigencomponents are given by the formula ( The meaning is consistent with that of ) The weights of the non-stationary characteristic components are given by the formula ( The meaning is consistent with that of ) The final predicted value is the result of weighted fusion of stationary and non-stationary feature components. With formula ( The meaning is consistent with that of ) This refers to the prediction results of the non-stationary feature sequences output by the DFCNformer model.

[0151] Step 4: Extract the receiver hardware differential code bias time series (DCB time series) from the BeiDou observation data, and use it as a rainfall prediction factor after preprocessing.

[0152] Specifically, such as Figure 3 As shown, firstly, the DCB time series of each receiver is extracted from the BeiDou observation data. Based on the GNSS dual-frequency ranging code and carrier phase observations, frequency-independent geometric distance, tropospheric delay, satellite clock bias, and receiver clock bias can be eliminated, resulting in the observation equation. Due to the significant noise in pseudorange observations, a carrier phase smoothing pseudorange method is used to extract DCB information, eliminating ionospheric delay. The sum of the satellite and receiver DCBs is taken as the observation value, as shown in the following formula:

[0153]

[0154] In the formula, For receiver differential code bias, For satellite differential code bias, For receiver and satellite In frequency and On The sum of observed values. In the absence of a constrained reference, the satellite and receiver... Since they cannot be separated, the sum of the DCBs of all satellites is usually set to zero, i.e., a "zero-mean" centroid reference is added. Then, through least squares adjustment, the DCBs of each receiver can be separated. Time series The specific formula is as follows:

[0155]

[0156] In the formula, The observation vector is the sum of the satellite and receiver DCB values. To design the matrix, For the satellite elevation angle weighting matrix, The normal equation matrix corresponding to the constraint datum. Let be the vector of parameters to be estimated. Then, use 3... Outliers were removed using conventional techniques (not detailed here), missing values ​​were imputed using linear interpolation (also a conventional technique), and high-frequency noise in the DCB time series was removed using moving average filtering (a conventional technique). Finally, the DCB sequence was subjected to singular spectrum analysis for feature separation, and stationary feature components were obtained. and Perform splicing, non-stationary feature components and The sets are then concatenated and divided into training, validation, and test sets.

[0157] Step 5: Obtain the vertical distribution data of atmospheric precipitable water from the reanalysis dataset, calculate the vertical distribution data of zenith wet delay and zenith tropospheric delay from the radiosonde data, and extract the vertical distribution characteristics of water vapor as a precipitation type determination factor.

[0158] Specifically, such as Figure 4 As shown, firstly, variables such as specific humidity (q), relative humidity (RH), temperature (t), and gravitational potential (z) at different pressure levels are downloaded from the ERA5 reanalysis dataset to obtain the vertical distribution data of PWV. Temperature, pressure, and specific humidity at each upper level are extracted from the University of Wyoming sounding database, and the ZWD and ZTD at each altitude are calculated through integration. Then, after unifying the elevation datum, a tomographic grid is defined. In the horizontal direction, bilinear interpolation is used to interpolate the two data points to the horizontal grid points of the target grid. In the vertical direction, exponential interpolation is used to resample the data to the target standard altitude level. The bilinear interpolation method assumes the interpolation point is... The four known points are The corresponding water vapor parameter values ​​are as follows: Furthermore, the interpolation point is located between these four known points, as shown in the following formula:

[0159]

[0160] In the formula, The value represents the water vapor parameter at the point to be interpolated.

[0161] Next, the vertical distribution characteristics of water vapor are calculated, including the vertical water vapor gradient. The vertical distribution center of water vapor The proportion of water vapor in the lower layers The proportion of water vapor in the middle layer The proportion of water vapor in the upper atmosphere Finally, all the calculated water vapor vertical distribution feature vectors were used as precipitation type determination factors.

[0162] Step 6: Construct and train a comprehensive rainfall prediction model. Use the final predicted values ​​of traditional prediction variables, the receiver hardware differential code deviation time series, and the water vapor vertical distribution characteristics as the prediction feature set, and then output the rainfall prediction result through the pre-trained comprehensive rainfall prediction model.

[0163] Specifically, the integrated rainfall prediction model combines classification and regression models. The dataset is divided into training, validation, and test sets. For the classification model, the vertical distribution feature vector of water vapor is used as the input variable. Based on a defined extreme rainfall threshold, historical extreme rainfall events are marked as "1," and a random forest algorithm is used to construct the rainfall discrimination model. The random forest consists of multiple decision trees, assuming a total number of trees... .

[0164] Decision trees, each tree For input features Make a prediction, either 0 or 1. The final prediction category is determined by the majority vote, using the following formula:

[0165]

[0166] In the formula, This represents the final predicted class of the random forest classification model, with a value of either 0 or 1. For majority voting functions, For the first The prediction function of a decision tree. , The total number of decision trees. The input feature vector is used. The Random Forest algorithm generates a training subset through bootstrap sampling and randomly selects a feature subset for splitting at each node to reduce inter-tree correlation and improve generalization performance. For the regression model, traditional predictor variables and DCB time series are used as input variables. The MixLinear model is used for prediction of stationary series, and the DFCNformer model is used for prediction of non-stationary series. The results are then weighted and fused to obtain the final result, which is then saved and output.

[0167] To verify the effectiveness of the proposed rainfall prediction method based on the fusion of two models (separation of stationary and non-stationary signal features), a comparative experiment was conducted. Long Short-Term Memory (LSTM) and Gate Recurrent Unit (GRU) networks were selected as baselines. These two methods are widely used models in the field of time series forecasting and are representative of existing GNSS rainfall prediction research. Observational data from the 244th to the 300th day of 2025 at an HKSC station in a certain region were used. The ZTD, ZWD, and PWV time series were obtained through step one. The comparative experiment used the exact same dataset, data preprocessing method, historical time window length (24 hours), and prediction step size (1 hour, 2 hours, and 3 hours) as the present invention. In both experiments, the samples were divided into a 60% training set, a 20% validation set, and a 20% test set in chronological order. LSTM and GRU were used to directly model and predict the original ZTD, ZWD, and PWV sequences, respectively. The evaluation metrics used are Mean Squared Error (MSE) and Mean Absolute Error (MAE). The smaller the values ​​of either metric, the higher the prediction accuracy. MSE is sensitive to large errors, while MAE reflects the average deviation of the predicted values. The two metrics are complementary and can comprehensively evaluate the model performance.

[0168] As shown in Tables 1 and 2, the MSE of this invention at prediction step sizes of 1 hour, 2 hours, and 3 hours are 0.102030, 0.213856, and 0.283528, respectively, and the MAE are 0.194530, 0.263615, and 0.317655, respectively. In terms of MSE, this invention improves upon the LSTM model by 72.09%, 63.65%, and 54.90%, respectively, and upon the GRU model by 72.32%, 64.40%, and 53.95%, respectively. In terms of MAE, this invention improves upon the LSTM model by 44.19%, 40.49%, and 34.27%, respectively, and upon the GRU model by 44.62%, 40.58%, and 33.94%, respectively.

[0169] Table 1. Comparison of prediction errors between the method of the present invention and LSTM and GRU.

[0170]

[0171] Table 2. Improvement rate of the method of the present invention compared with LSTM and GRU

[0172]

[0173] The results show that the present invention significantly outperforms LSTM and GRU models at all prediction step lengths, and the improvement in MSE is greater than that in MAE, indicating a more significant advantage in capturing extreme biases. This further demonstrates that by separating stationary and non-stationary features from the GNSS atmospheric parameter time series and employing targeted prediction models for each, the accuracy and stability of rainfall prediction can be effectively improved. In summary, through comparative experiments with two representative prediction methods, LSTM and GRU, the present invention achieves optimal prediction accuracy at prediction step lengths of 1h, 2h, and 3h, fully demonstrating the effectiveness and superiority of the technical solution of the present invention.

[0174] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely preferred examples and are not intended to limit the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A rainfall prediction method based on the fusion of dual models separating stationary and non-stationary signal features, characterized in that, Includes the following steps: Step 1: Acquire and preprocess rainfall forecasting-related data, which includes BeiDou raw observation data, historical radiosonde data, and reanalysis datasets; Step 2: Separate stationary and non-stationary feature sequences from the preprocessed BeiDou raw observation data using singular spectrum analysis. Step 3: Construct a first prediction model to predict the stationary feature sequence and obtain the stationary feature prediction result; construct a second prediction model to predict the non-stationary feature sequence and obtain the non-stationary feature prediction result; weight and fuse the stationary feature prediction result and the non-stationary feature prediction result to obtain the final predicted value of the traditional predictor variable. Step 4: Extract the receiver hardware differential code deviation time series from the BeiDou observation data, and use it as a rainfall prediction factor after preprocessing. Step 5: Obtain the vertical distribution data of atmospheric precipitable water from the reanalysis dataset, calculate the vertical distribution data of zenith wet delay and zenith tropospheric delay from the radiosonde data, and extract the vertical distribution characteristics of water vapor as a precipitation type determination factor. Step 6: Construct and train a comprehensive rainfall prediction model. Use the final predicted values ​​of traditional prediction variables, the receiver hardware differential code deviation time series, and the water vapor vertical distribution characteristics as the prediction feature set, and then output the rainfall prediction result through the pre-trained comprehensive rainfall prediction model.

2. The method according to claim 1, characterized in that, The specific steps for acquiring and preprocessing rainfall forecast data include: The original BeiDou observation data was processed using precise point positioning technology, the cutoff elevation angle was set, the ionospheric delay effect was eliminated by using an ionospheric-free combination model, and the zenith tropospheric delay was estimated by using broadcast ephemeris files and precise orbit and clock difference files through Kalman filtering. Based on the air pressure, latitude, and altitude at the receiver site, the zenith static delay is calculated using the Sastamonin model. The zenith wet delay is obtained by subtracting the zenith static delay from the zenith tropospheric delay. Precipitation is calculated based on the zenith wet delay and the conversion factor. The zenith tropospheric delay, zenith wet delay, and precipitation constitute traditional rainfall prediction variables. Vertical profile data of temperature, humidity and air pressure are extracted from the sounding database, and atmospheric precipitable water and pressure layer data are obtained from the reanalysis dataset to obtain water vapor vertical distribution data.

3. The method according to claim 1, characterized in that, In step 2, the preprocessed zenith tropospheric delay time series, zenith wet delay time series, and precipitable water time series are decomposed into time series to separate stationary and non-stationary characteristic series.

4. The method according to claim 3, characterized in that, In step 2, the specific steps for separating stationary and non-stationary feature sequences include: Singular spectral analysis is used for feature separation. The window length is set, and the time series is mapped to a trajectory matrix by an embedding operator. Singular value decomposition is performed on the trajectory matrix to obtain multiple components; Components associated with steady trends and periodicity are grouped into one group, while components associated with noise and non-stationary abrupt changes are grouped into another group. The two components are reconstructed into one-dimensional time series using the diagonal averaging method, yielding stationary and non-stationary characteristic sequences.

5. The method according to claim 1, characterized in that, The specific steps for constructing the first prediction model to predict stationary feature sequences and obtaining stationary feature prediction results include: The first prediction model is a hybrid linear prediction model, which includes a time-domain segmented processing path and a frequency-domain adaptive low-rank filtering path. The stationary feature sequence is input into the hybrid linear prediction model. Local waveform features are extracted and cross-segment dependencies are captured through the time-domain segmentation processing path. The downsampled data sequence is then subjected to Fourier transform and low-rank decomposition adaptive spectral filtering through the frequency-domain adaptive low-rank filtering path. The outputs of the two paths are upsampled and then fused by addition to obtain the prediction results of the stationary feature components.

6. The method according to claim 5, characterized in that, Extracting local waveform features and capturing cross-segment dependencies through time-domain segmentation specifically includes: The input sequence is downsampled and divided into multiple non-overlapping time segments; Perform intra-segment linear projection on each time segment to extract local waveform features; After concatenating all intra-segment features, perform inter-segment linear projection to capture cross-segment dependencies.

7. The method according to claim 5, characterized in that, The adaptive low-rank filtering path in the frequency domain performs Fourier transform and low-rank decomposition adaptive spectral filtering on the downsampled data sequence, specifically including: Apply the Fast Fourier Transform to the downsampled data sequence to obtain the Fourier Transform result; Adaptive spectral filtering is achieved by employing low-rank decomposition technology. The spectrum of each frequency band is mapped to a low-dimensional latent space and reconstructed to obtain the frequency domain representation after low-rank adaptive spectral filtering.

8. The method according to claim 1, characterized in that, The specific steps for constructing a second prediction model to predict the non-stationary feature sequence and obtain the non-stationary feature prediction results include: The second prediction model is a destationary Fourier transform and coefficient network prediction model; The non-stationary feature sequences are normalized and stabilized to ensure that the statistical properties of the input data remain consistent. The processed sequences were decomposed into seasonal and trend components; For the seasonal components, a de-stationary Fourier attention mechanism is used, and an encoder-decoder architecture is employed for prediction to obtain the prediction results of the seasonal components. The trend components are predicted using a multilayer perceptron and a dual-coefficient network to obtain the prediction results. The prediction results of seasonal components are added together with the prediction results of trend components, and then inverse stationarization is performed to obtain the prediction results of non-stationary characteristic components.

9. The method according to claim 8, characterized in that, The processed sequence is decomposed into seasonal and trend components, specifically including: By using moving average filters with different window sizes to capture various trend patterns, and by using adaptive weights to merge all trend patterns, the final trend component is obtained. The final seasonal component is obtained by subtracting the final trend component from the original processed sequence.

10. The method according to claim 8, characterized in that, For the seasonal components, the specific steps for prediction using a de-stationary Fourier attention mechanism and an encoder-decoder architecture include: The seasonal components are processed through a linear layer to obtain a query matrix, a key matrix, and a value matrix. Calculate the non-stationary factor, which includes the standard deviation of the input sequence and the mean value of the query matrix over the time dimension; The query matrix, key matrix, and value matrix are transformed into the frequency domain using Fourier transform; An attention mechanism is applied in the frequency domain, and the non-stationary factor is used to correct the attention calculation results; The frequency domain results are converted back to the time domain by inverse Fourier transform, and then passed to the multilayer encoder and multilayer decoder for processing to obtain the prediction results of the seasonal components.

11. The method according to claim 8, characterized in that, The specific steps for predicting the trend components using a multilayer perceptron and a dual-coefficient network include: The trend component sequence is passed through the input coefficient network to obtain the horizontal parameters and fluctuation parameters of the input data, and then the output is normalized. The normalized output is input into the multilayer perceptron to complete model training; After the multilayer perceptron outputs the prediction result, it predicts future values ​​through the output parameter network and outputs the normalized prediction result. The normalized prediction results are then subjected to inverse normalization to obtain the prediction results of the trend components.

12. The method according to claim 1, characterized in that, The specific steps for weighted fusion of the stationary feature prediction results and the non-stationary feature prediction results to obtain the final predicted value of the traditional predictor variable include: The weights of stationary and non-stationary feature components are determined by using Bayesian optimization methods to minimize the prediction error of the validation set. Based on the weights of the stationary feature components and the non-stationary feature components, the prediction results of the stationary feature and the non-stationary feature are weighted and summed to obtain the final predicted value of the traditional predictor variable, wherein the sum of the weights of the stationary feature components and the non-stationary feature components is 1.

13. The method according to claim 2, characterized in that, The specific steps for extracting the receiver hardware differential code offset time series from BeiDou observation data and preprocessing it as a rainfall prediction factor include: Based on the dual-frequency ranging code and carrier phase observations, the differential code deviation information is extracted using the carrier phase smoothing pseudorange method to eliminate ionospheric delay and obtain the observed value of the sum of satellite differential code deviation and receiver differential code deviation. Add a zero-mean centroid reference and separate the differential code deviation time series of each receiver by least squares adjustment; Outliers were removed using statistical methods, missing values ​​were imputed using linear interpolation, and high-frequency noise was removed using moving average filtering to obtain the preprocessed receiver hardware differential code deviation time series. The receiver hardware differential code bias time series is decomposed into stationary and non-stationary feature components, which are then concatenated with the stationary and non-stationary feature sequences of the traditional rainfall prediction variables, respectively, and used as input to the prediction model.

14. The method according to claim 1, characterized in that, The specific steps for obtaining vertical distribution data of atmospheric precipitable water from the reanalysis dataset, calculating vertical distribution data of zenith wet delay and zenith tropospheric delay from radiosonde data, and extracting water vapor vertical distribution characteristics as precipitation type determination factors include: Download the specific humidity, relative humidity, temperature and gravitational potential variables of different pressure layers from the reanalysis dataset to obtain the vertical distribution data of atmospheric precipitable water. Temperature, air pressure and specific humidity of each altitude layer are extracted from the sounding database, and the zenith wet delay and zenith tropospheric delay of each altitude layer are calculated by integration. After unifying the elevation datum, a tomographic grid is defined. In the horizontal direction, bilinear interpolation is used to interpolate the data to the horizontal grid points of the target grid, and in the vertical direction, exponential interpolation is used to resample the data to the target standard height layer. The vertical distribution characteristics of water vapor are calculated based on the interpolated data. These characteristics include the vertical gradient of water vapor, the height center of the vertical distribution of water vapor, the proportion of water vapor in the lower layer, the proportion of water vapor in the middle layer, and the proportion of water vapor in the upper layer.

15. The method according to claim 1, characterized in that, The specific steps for constructing a comprehensive rainfall prediction model include: The samples are divided into training set, validation set and test set in chronological order, and the historical time window length, sampling interval and prediction step size are set. Construct a comprehensive rainfall prediction model that combines classification and regression models; For the classification model, the vertical distribution characteristics of water vapor are used as input variables. Historical extreme rainfall events are labeled as specified categories according to the set extreme rainfall threshold. A random forest algorithm is used to construct a rainfall discrimination model, and the final predicted category is determined by majority voting. For the regression model, the traditional rainfall prediction variables and the receiver hardware differential code deviation time series are used as input variables. The first prediction model is used to predict stationary series, and the second prediction model is used to predict non-stationary series. The prediction results are then weighted and fused to obtain the final rainfall prediction result. The discrimination result of the classification model and the prediction result of the regression model are used together as the output of the integrated rainfall prediction model.

16. The method according to claim 15, characterized in that, The random forest algorithm generates a training subset by sampling with replacement and randomly selects a feature subset at each node for splitting, thereby reducing the correlation between trees and improving generalization performance.

Images