A method for identifying, spatially classifying, and analyzing driving factors of rainstorm events

Abstract

This invention discloses a method for identifying, spatially classifying, and analyzing driving factors of rainstorm events. The method includes: acquiring and processing meteorological observation and reanalysis data to obtain standardized precipitation time series and synchronous meteorological element fields; employing a physically constrained dynamic time window mechanism to adaptively identify and classify rainstorm events, generating a physically consistent rainstorm event sequence set; using a Bayesian fusion framework to integrate multiple statistical test evidences and physical prior information based on the climate oscillation index to detect abrupt change characteristics of the sequence set and generate a high-confidence abrupt change point set; and finally, integrating and analyzing the data to generate a report. This invention improves the physical accuracy of rainstorm event identification and the reliability of abrupt change detection.

Description

Technical Field

[0001] This invention belongs to the field of meteorological and hydrological data processing, and in particular, it is a method for identifying, spatially classifying, and analyzing driving factors of rainstorm events. Background Technology

[0002] Against the backdrop of global climate change, the frequency, intensity, and spatiotemporal distribution of rainstorm events exhibit significant non-stationarity. In-depth analysis of the evolution patterns of rainstorm events, accurate identification of abrupt changes in their long-term trends, and understanding of the underlying physical driving mechanisms are of crucial scientific value and practical significance for improving meteorological forecasting and early warning capabilities, formulating scientific flood control and disaster reduction strategies, and conducting sustainable water resource management.

[0003] Currently, research and analysis of rainstorm events have made some progress. Existing techniques typically rely on daily precipitation observation data collected from ground meteorological stations. First, a fixed precipitation threshold (e.g., daily precipitation greater than or equal to 50 mm) is set to identify single-day rainstorms. Based on this, a fixed time interval without precipitation (e.g., 24 or 48 consecutive hours without effective precipitation) is defined as the standard for event segmentation to divide independent rainstorm processes and statistically analyze their annual characteristics such as frequency and total precipitation. For these constructed rainstorm characteristic time series, researchers often use non-parametric statistical methods such as the Mann-Kendall (MK) test to detect the presence of long-term trends and abrupt changes in the series, and combine this with traditional tools such as wavelet analysis to identify their periodic patterns.

[0004] However, existing technologies suffer from a deep-seated technical problem when dealing with phenomena like rainstorms, which have complex physical causes: a disconnect between statistical definitions and physical reality. This is mainly manifested in the lack of physical consistency in event segmentation, leading to insufficient reliability of mutation detection results. Summary of the Invention

[0005] The purpose of this invention is to provide a method for identifying, spatially classifying, and analyzing driving factors of rainstorm events, in order to solve the aforementioned problems existing in the prior art.

[0006] Technical solution: A method for identifying, spatially classifying, and analyzing driving factors of rainstorm events, including:

[0007] Acquire and process meteorological observation and reanalysis data to obtain standardized precipitation time series and synchronous meteorological element fields;

[0008] Based on standardized precipitation time series and synchronous meteorological element fields, a physical constraint-based dynamic time window mechanism is adopted to adaptively identify rainstorm events and generate a rainstorm event sequence set.

[0009] Based on the rainstorm event sequence set and externally acquired climate oscillation index data, a Bayesian fusion framework is used to integrate statistical test evidence and physical prior information to detect abrupt change characteristics of the rainstorm event sequence set and generate a high-confidence abrupt change point set.

[0010] By integrating rainstorm event sequence sets, high-confidence abrupt change point sets, and standardized precipitation time series, multi-dimensional feature analysis is performed to generate a comprehensive analysis report on rainstorm characteristics.

[0011] Beneficial effects: This invention improves the physical accuracy of rainstorm event identification and the reliability of mutation detection. Attached Figure Description

[0012] Figure 1 A flowchart illustrating the steps of a method for identifying, spatially classifying, and analyzing driving factors of rainstorm events, provided in an embodiment of this application.

[0013] Figure 2 A flowchart illustrating the steps for generating a rainstorm event sequence set as provided in this application embodiment.

[0014] Figure 3 A flowchart illustrating the steps for forming a weather system transition time set provided in this application embodiment.

[0015] Figure 4 A flowchart illustrating the steps for generating a set of conversion points after physical verification, as provided in an embodiment of this application. Detailed Implementation

[0016] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0017] It should be noted that the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion, for example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0018] The study found that the currently widely used fixed no-precipitation intervals as the sole basis for event segmentation are arbitrary in terms of physical mechanisms. They fail to distinguish between two fundamentally different situations: first, two temporally segmented precipitation events caused by the same slowly moving or intermittently developing weather system may be incorrectly separated into two independent events; second, two weather systems with completely different physical origins may be incorrectly merged into a single event if they happen to pass through close together. This bias in event definition undermines the physical consistency of the source data, causing the constructed rainstorm event sequence set to fail to accurately reflect the formation and dissipation processes of weather systems. Furthermore, the reliability of subsequent abrupt change tests (such as the MK test) decreases when applied to sequences with such physical consistency biases and strong autocorrelation. The test method may identify statistically significant abrupt changes, but these abrupt changes may not physically correspond to any real shifts in the large-scale climate background, but are merely statistical artifacts introduced by incorrect event segmentation, leading to a high misclassification rate for years of climate abrupt changes.

[0019] An exemplary scheme for analyzing the characteristics of rainstorm events includes the following steps: collecting multi-year daily precipitation observation data from multiple meteorological stations in the target watershed and performing quality control; identifying rainstorm events based on a fixed threshold method (e.g., daily precipitation greater than or equal to 10 mm) and extracting characteristic indicators such as the number of rainstorm days, total precipitation, and number of rainstorm events; analyzing the trend and abrupt changes of the sequence using the Theil-Sen and MK tests, and extracting periodic features using wavelet packet decomposition; classifying rainstorm events and analyzing their relationship with atmospheric circulation factors.

[0020] Specifically, multi-year daily precipitation observation data for the target watershed are collected, and missing and outlier values ​​are checked and corrected to ensure the integrity and consistency of the time series, forming a complete precipitation data sequence. A daily precipitation database containing multiple meteorological stations is established. Heavy rain events in the daily precipitation sequence are identified using a fixed threshold method. The determination rule is: when the daily precipitation is greater than or equal to a set threshold (e.g., 10 mm / d), that day is a heavy rain day; if the threshold condition is met for several consecutive days, it constitutes a heavy rain event. For each meteorological station, the following rainstorm characteristic indicators are extracted annually: number of rainstorm days, total rainstorm amount, number of rainstorm events, and maximum annual rainstorm amount. This forms a multi-dimensional attribute feature sequence for rainstorms, and a rainstorm feature sequence dataset is constructed. The specific process is as follows: If the precipitation on a certain day exceeds the threshold of 10 mm, it is counted as one rainstorm day. The total number of rainstorm days is accumulated throughout the year (statistically on an annual scale). The cumulative value of the precipitation on all rainstorm days is used as the total rainstorm amount for that year (statistically on an annual scale). Multiple consecutive days exceeding the threshold are judged as the same rainstorm event. When the precipitation falls below the threshold, the rainstorm event is considered to have ended. The characteristics of the event are recorded, and the identification of the next event begins. The total number of rainstorm events accumulated annually is the number of rainstorm events for the corresponding year. The sum of the precipitation on all rainstorm days during a rainstorm is the amount of rainstorm event. The maximum value of multiple rainstorm events each year is the maximum annual rainstorm amount. Non-stationarity detection and multi-scale analysis were performed on the multidimensional characteristic sequences of rainstorms, including: identifying the long-term evolution direction of rainstorm characteristics using Theil-Sen estimation coupled with the M-K test, and based on the trend-preserving pre-white sequence x*. t The standardized statistic Z is used to determine trend and significance; the positive order statistic UF is calculated using the M-K mutation test. k With reverse order statistic UB k The year corresponding to the intersection of the two is taken as a possible abrupt change point; wavelet packet decomposition is used to analyze the periodicity of the rainstorm characteristics to determine the main period T*. Further, specifically including: Let the original time series be: {x} t}, t = 1, ..., T. The Theil–Sen method is used to obtain the median slope between every two data points as the slope estimate: β* = median{(x j -x i ) / (ji), 1≤i<j≤T}; the intercept estimate is taken as the median residual: α*= median{x t -β*t};Detrend, obtain the residual sequence: r t = x t -(α*+β*t); using the Yule–Walker formula in the residual r t Estimate the autoregressive coefficient Φ*: Φ* = ∑ t=2 T r t rt-1 / ∑ t=1 T-1 r t 2 Construct a trend-preserving pre-white sequence x* t :x* t = x t -Φ* x t-1 , t=2, …,T; where the first term is treated as x*1= x1. At this point, from x t To x* t The transformation weakens the autocorrelation component, while the Φ* estimated on the detrended residuals preserves the linear trend of the original time series. In x* t Calculate the M-K statistic S and the standardized statistic Z: S = ∑ i=2 n ∑ j=1 i-1 sgn(x* i - x* j ), sgn(x* i - x* j )=-1,x* i - x* j <0; 0, x* i - x* j =0; 1, x* i - x* j >0; Z=sqrt(18)(S-1) / sqrt(n(n-1)(2n+5)), S>0; 0, S=0; sqrt(18)(S+1) / sqrt(n(n-1)(2n+5)), S<0; where n is the sample size; a Z value between 0 and 1.96 indicates an increasing trend but the trend is not significant; a value between 1.96 and 2.64 indicates that at a significance level α=0.05, the 95% confidence level significance test is passed, indicating a significant increasing trend; if the value is greater than 2.64, it indicates that the significance test is greater than 99%, indicating a highly significant increasing trend. Based on x* t The M-K statistic S is calculated above, and the standardized statistic UF of the series is calculated. k With the inverse sequence statistic UB k This involves using an improved M-K mutation test to identify mutation points in the characteristic sequences of rainstorms. For sequence x* t Perform maximum resolution wavelet packet decomposition (WP) to determine the main period. Specifically, this includes: for sequence x... t Perform wavelet packet decomposition down to layer J to obtain the coefficient sequence c for each node (J, n). J，n [k]. Calculate the subband energy and normalized energy distribution: E J，n =∑ kc J，n [k] 2 E tot =∑ n=0 2J-1 E J，n ;p J，n = E J，n / E tot Among them, E J，n The subband energy, i.e., the sum of squares of the coefficients; E tot p represents the total energy. J，n This represents the energy percentage of the subband. Calculate the energy spectral entropy: H = -∑ n=0 2J-1 p J，n ln(p J，n ); H norm =H / ln(2 J H ∈ [0,1]; where H is the Shannon energy spectrum entropy; H norm The normalizable entropy is a smaller value indicating more concentrated energy, i.e., a more distinct principal period, while a larger value indicates more dispersed energy. Extracting the principal period T*: n* = arg max n p J，n T*=T J，n* p*= p J，n* Where n* is the index of the largest energy sub-band, which is the main energy sub-band number, i.e., the energy percentage of all sub-bands p. J，n In the process, find the sub-band number with the largest energy percentage; T* is the main period, corresponding to the center period of the main energy sub-band, and is determined by the center frequency f of that sub-band. J，n*The converted value represents the most significant periodic change cycle in the time series; p* represents the energy proportion corresponding to the main cycle, indicating the relative importance of this cycle in the total energy. The higher the value, the more dominant the cycle. By analyzing the energy distribution and main cycle of different sub-bands, the dominant oscillation scale and secondary periodic components of the rainstorm characteristic sequence can be revealed. A graded threshold method is used to classify rainstorm events according to the total amount and duration of rainfall: in terms of magnitude, they are divided into rainstorm, heavy rainstorm, and extremely heavy rainstorm; in terms of duration, they are divided into short-duration rainstorm and long-duration rainstorm. By statistically analyzing the spatiotemporal distribution patterns of rainstorms of different grades, the spatial heterogeneity within the watershed is characterized, and typical events are selected for comparative analysis. Furthermore, the graded thresholds for magnitude are set as: 50 mm; 250 mm, dividing them into rainstorm (50–99.9 mm), heavy rainstorm (100–249.9 mm), and extremely heavy rainstorm (≥250 mm). By classifying rainstorms into heavy, torrential, and extremely heavy intensities, the spatial distribution characteristics and regional differences of rainstorms of varying intensities within the watershed can be revealed. A duration threshold of 3 days was set, dividing rainstorms into short-duration (1–2 days) and long-duration (≥3 days) rainstorms. Distinguishing between short-duration and long-duration rainstorms reflects the differences in the duration of rainstorm events, providing a basis for identifying regionally severe convective short-duration rainstorms and large-scale long-duration rainstorms. Typical representative stations in the watershed were selected to analyze the correlation between rainstorm characteristic sequences and major atmospheric circulation factors (including the North Atlantic Oscillation (NAO), Arctic Oscillation (AO), East Asian Summer Monsoon (EASM), and Antarctic Oscillation (AAO). The linear correlation was determined using the Pearson correlation coefficient r; the correlation was further analyzed using cross-wavelet transform W. xy Calculate the cross wavelet power |W (a, b) xy (a, b)∣ reveals the significant coupling relationship between rainstorms and climate factors at different time scales, and determines their leading-lag effect through phase difference. Furthermore, Pearson correlation analysis is used to reveal the linear correlation between rainstorm characteristic sequences and atmospheric circulation factors. The correlation coefficient r is calculated as: r = ∑ i=1 n (x i -x')(y i -y') / sqrt(∑ i=1 n (x i -x') 2 )sqrt(∑ i=1 n (y i -y') 2 ); where x i For the characteristic sequence of rainstorms, y iLet x' and y' be the climate factor sequence, x' and y' be the means, and n be the sample length. Cross-wavelet correlation analysis is used to identify the coupling relationship between the rainstorm characteristic sequence and the atmospheric factor sequence in the time-frequency domain. Let the rainstorm sequence be x(t) and the atmospheric factor sequence be y(t), and their wavelet transforms be W... x (a, b) and W y (a, b), then the cross wavelet spectrum is defined as: W xy (a, b) = W x (a, b) W y *(a, b), where * denotes complex conjugate. The corresponding cross wavelet power is: |W xy (a, b)∣=∣W x (a, b) W y *(a, b)∣; By using the significance test of cross-wavelet power, we can identify the significant coupling relationship between rainstorm characteristics and atmospheric circulation factors at different scales and times; at the same time, we can reveal the lead-lag relationship between the two through phase difference analysis.

[0021] In this embodiment, a fixed threshold method is used to identify and classify rainstorm events, for example, classifying multi-day events exceeding a threshold as the same rainstorm event. This method, based on fixed time intervals (e.g., 24 hours without rain) or fixed precipitation thresholds, is physically crude and ignores the inherent physical continuity of the weather system between two precipitation events. In other words, two precipitation events separated in time may still belong to the same large, slowly moving, or intermittently developing weather system, while the fixed threshold method incorrectly separates them into two independent rainstorm events. Conversely, two short-duration heavy precipitation events with completely different physical causes, if the time interval is close, may also be misclassified as the same event. This bias in classification affects the statistical accuracy of subsequent rainstorm frequency, intensity, and other characteristics. Furthermore, the MK test is used to determine abrupt change points. Single statistical tests are sensitive to the autocorrelation of time series and are prone to producing false positive abrupt change points in periods without real physical background changes, leading to misjudgments of years of climate change. Analyzing precipitation data on a daily scale makes it impossible to capture the fine structure within a rainstorm event, such as key features like hourly peak rainfall intensity and the time of peak rainfall, resulting in an incomplete characterization of the extreme nature of rainstorms.

[0022] Therefore, such as Figure 1 As shown, a method for identifying, spatially classifying, and analyzing driving factors of rainstorm events is proposed, including:

[0023] Meteorological observation and reanalysis data are acquired and processed to obtain standardized precipitation time series and synchronous meteorological element fields.

[0024] In other words, by acquiring meteorological observation and reanalysis data and performing preprocessing, standardized precipitation time series and synchronous meteorological element fields can be obtained.

[0025] Specifically, meteorological observation data can include, but is not limited to, hourly or daily precipitation observation sequences recorded by surface meteorological stations. Reanalysis data can be a spatiotemporally continuous and physically consistent gridded meteorological dataset generated by combining observational data from multiple sources (such as surface, radiosonde, and satellite) with numerical weather prediction models, such as the ERA5 reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF). Standardized precipitation time series refer to precipitation sequences that have undergone quality control, interpolation, autocorrelation removal, and standardization (e.g., Z-score standardization) to facilitate subsequent statistical analysis. Synchronous meteorological element fields refer to multidimensional meteorological physical quantity fields that completely correspond to the precipitation time series in time and space, such as elements including potential temperature, specific humidity, and wind field on specific isobaric surfaces.

[0026] Based on standardized precipitation time series and synchronous meteorological element fields, a physical constraint-based dynamic time window mechanism is adopted to adaptively identify rainstorm events and generate a rainstorm event sequence set.

[0027] In other words, based on standardized precipitation time series and synchronous meteorological element fields, by coupling integrated meteorological similarity calculation, change point detection algorithm and water vapor flux continuity verification, a dynamic time window mechanism is used to achieve adaptive identification of rainstorm events and output a rainstorm event sequence set.

[0028] In this embodiment, the traditional method of dividing rainstorm events using fixed time intervals (such as 24 hours) or fixed precipitation thresholds is abandoned. Instead, physical information (such as water vapor flux, thermodynamic parameters, etc.) in a synchronous meteorological element field is used to determine whether adjacent precipitation processes belong to the same weather system, and the length of the precipitation-free time window used to divide the events is dynamically adjusted based on information such as precipitation intensity and system movement speed. This results in higher consistency and completeness in the physical causes of each identified rainstorm event. For example, the rainstorm event sequence set contains a collection of multiple identified independent rainstorm events, each with a clearly defined start and end time and other attributes.

[0029] Based on the rainstorm event sequence set and externally acquired climate oscillation index data, a Bayesian fusion framework is used to integrate statistical test evidence and physical prior information to detect abrupt change characteristics of the rainstorm event sequence set and generate a high-confidence abrupt change point set.

[0030] In other words, based on the rainstorm event sequence set and the climate oscillation index data, the detection evidence of multiple mutation testing methods and the physical prior constructed based on the climate oscillation index data are integrated through the Bayesian posterior probability framework to detect the mutation characteristics of the rainstorm sequence and generate a high-confidence mutation point set and a mutation intensity assessment matrix.

[0031] Specifically, heavy rainfall event sequences can be processed into interannual scale characteristic sequences, such as annual heavy rainfall frequency or annual total heavy rainfall. Climate oscillation index data refer to indicators characterizing the interannual or interdecadal variations of large-scale climate systems, such as the ENSO index (e.g., Nino 3.4), the Pacific Decadal Oscillation (PDO) index, and the Atlantic Multidecadal Oscillation (AMO) index. The Bayesian fusion framework can combine evidence from multiple sources (i.e., results from various traditional abrupt change testing methods, such as Mann-Kendall and Petittt) with prior knowledge based on physical mechanisms (i.e., whether large-scale climate background changes support the occurrence of abrupt changes) to calculate the posterior probability of each candidate abrupt change point. For example, a high-confidence abrupt change point set refers to the set of years whose posterior probabilities exceed a certain high threshold (e.g., 0.8); these abrupt changes are more reliable because they are supported by both statistical evidence and physical priors.

[0032] By integrating rainstorm event sequence sets, high-confidence abrupt change point sets, and standardized precipitation time series, multi-dimensional feature analysis is performed to generate a comprehensive analysis report on rainstorm characteristics.

[0033] In this embodiment, multi-dimensional feature analysis may include, but is not limited to: analyzing the long-term trends of rainfall frequency and intensity (e.g., using Sen's slope estimation); analyzing its multi-scale periodic characteristics (e.g., using wavelet analysis); analyzing its feature differences before and after different abrupt change points; and analyzing its spatial distribution pattern (e.g., generating heat maps through spatial interpolation). For example, a comprehensive rainfall feature analysis report may be a document containing a series of visualizations, statistical parameters, and conclusive descriptions, providing a scientific basis for climate change research and flood prevention and disaster reduction decision-making.

[0034] This embodiment improves the physical consistency of rainstorm event identification and the reliability of mutation detection through high-quality data processing, physical mechanism-based event recognition, multi-evidence fusion mutation detection, and comprehensive feature analysis.

[0035] In an optional embodiment, a standardized precipitation time series and synchronous meteorological element field are obtained, specifically including:

[0036] Hourly precipitation observation sequences (e.g., P) in meteorological observation and reanalysis data raw (t)) The quality control process based on multi-criteria identification and adaptive interpolation is applied to generate the post-quality control precipitation sequence (P).qc (t)).

[0037] Preferably, generating a post-quality control precipitation sequence includes:

[0038] By comprehensively applying four criteria—physical upper limit, statistical deviation, temporal consistency, and spatial consistency—the hourly precipitation observation sequence is examined to generate an outlier marker sequence.

[0039] In other words, multiple criteria are applied to identify and test outliers, generating a sequence of outlier labels.

[0040] Specifically, multiple physical and statistical constraints are used to identify suspicious values ​​in the raw data. Examples include: physical upper limit test: marking data exceeding local or physically possible precipitation extremes (e.g., 300 mm / h) as outliers; statistical deviation test: calculating the mean μ over a relatively long (e.g., 72-hour) sliding window. w and standard deviation σ w and mark deviations from μ w ±3σ w Data at (or other selected multiples) are statistically abnormal; temporal consistency test: check the instantaneous rate of change of the precipitation series. If the precipitation jumps more than a large threshold (e.g., |P|) between adjacent time points (e.g., within 1 hour), the result is considered statistically abnormal. raw (t) - P raw If the precipitation at a station is greater than 100 mm / h and the jump is not sustained (e.g., less than 2 hours), it is marked as an abrupt change anomaly. Spatial consistency test: The precipitation value at this station is compared with the average value of the same period at surrounding (e.g., 5) neighboring stations. If the difference exceeds a certain threshold (e.g., 50 mm / h), it is marked as a spatial anomaly. By combining one or more of the above criteria, an outlier labeling sequence Flag is generated. anomaly (t).

[0041] Based on outlier labeling sequence analysis, the temporal distribution characteristics of missing data are classified into different missing data patterns.

[0042] In other words, based on the null values ​​and marked outliers in the raw data, we can analyze the duration and patterns of consecutive missing data. For example, they can be categorized as: short-term missing data (e.g., <3 hours), medium-term missing data (e.g., 3-12 hours), long-term missing data (>12 hours), and systematic missing data (e.g., occurring at specific times each day).

[0043] For any missing data pattern, a matching interpolation strategy is selected from the preset interpolation strategy library to fill in the data and obtain the quality-controlled precipitation sequence.

[0044] In this embodiment, adaptive interpolation is implemented, meaning that different and most suitable mathematical methods are used for different missing data patterns. Preferably, for short-term missing data, cubic spline interpolation can be used to ensure the smoothness of the interpolation sequence and the continuity of its derivatives. For moderate missing data, linear interpolation can be used to obtain preliminary values, which are then corrected using the hourly precipitation increments from adjacent stations during the same period. For long-term missing data, a weighted average of neighboring stations can be used, where the weight w... i It can be determined by distance d i and correlation coefficient r i Joint decision (e.g., w) i =r i / (d i 2 For systematic missing data, in addition to the spatial interpolation mentioned above, historical climatological data from the same period can be used as a reference for correction. Through the above interpolation strategy, the preliminarily filled interpolated sequence P is obtained. filled (t). For the interpolated sequence P filled (t) Perform a rationality check and output the quality-controlled precipitation sequence P. qc (t). Necessary verification to ensure that the interpolation results do not introduce new errors. For example, reasonableness checks may include: physical reasonableness (e.g., the interpolation result must be non-negative and not exceed a specific multiple of the region's historical extreme values, such as 1.5 times); temporal continuity (the connection between the interpolated segment and the preceding and following measured segments should be smooth, and the first derivative should not have extreme jumps); and statistical consistency (the change in monthly or annual cumulative precipitation before and after interpolation should be controlled within a certain proportion, such as 20%).

[0045] For the quality control precipitation sequence P qc (t) Perform segmented pre-whitening processing (to remove autocorrelation), and incorporate a trend information preservation and verification mechanism into the segmented pre-whitening processing. Then, standardize the quality-controlled precipitation series after segmented pre-whitening processing (e.g., Z-score standardization) to obtain the standardized precipitation time series P. std (t). Among them, the pre-whitened sequence (P) may have undergone trend correction. w_seg (t) or P wtrend (t) is used as input for subsequent trend analysis (i.e., the pre-whitening precipitation sequence P). w (t)).

[0046] In this embodiment, the interference of the temporal correlation of the precipitation sequence itself on subsequent trend and abrupt change detection (such as the MK test) is eliminated, while ensuring that the true climate change trend signal is not erroneously rejected. Specifically: for the quality-controlled precipitation sequence P qc(t) Perform the augmented Dickey-Fuller test and the KPSS test to determine the stationarity of the series and the possible differencing order d. If the series is non-stationary, perform differencing of order d first. Calculate the partial autocorrelation function (PACF) of the differrated series and identify the maximum lag order p with PACF truncation. max From 1 to p max Fit autoregressive (AR) models of different orders within a given range, calculate the Bayesian information criterion (BIC) for each model, and select the order that minimizes the BIC as the optimal AR order p for the sequence. opt .

[0047] In a preferred implementation, the post-quality control precipitation sequence undergoes segmented pre-whitening processing to remove autocorrelation, and a trend information preservation verification mechanism is incorporated into the segmented pre-whitening processing, including:

[0048] The quality-controlled precipitation sequence is divided into predetermined segments according to age. For any segment, an autoregressive model is independently fitted to construct the pre-whitening transformation of that segment, generating a segmented pre-whitening sequence. The long-term trends of the segmented pre-whitening sequence and the quality-controlled precipitation sequence are calculated separately, and the differences between the two are compared. If the difference exceeds a preset threshold, a trend-preserving term is introduced into the pre-whitening transformation to correct the segmented pre-whitening sequence.

[0049] Specifically, the quality-controlled precipitation sequence P qc (t) is divided into multiple sub-segments according to time (e.g., every 10 years). For any given sub-segment, its AR(p) is independently fitted. opt The model estimates the autoregressive coefficients Φ1, Φ2, ..., Φ3 specific to this time period. p The Z-score standardization of this segment yields the standardized precipitation time series P. std (t), construct the pre-whitening transformation for this segment: P w_seg (t)=P std (t)-∑ i=1 p Φ i ×P std (ti), generate segmented pre-whitening sequence P w_seg (t). This segmented processing allows the pre-whitening model to adapt to different interdecadal climate characteristics. The segmented pre-whitening sequences P are calculated separately. w_seg The long-term trend of (t) (e.g., using Sen's slope estimation to obtain β) whitened ) and (unpre-bleached) quality control precipitation sequence P qc The long-term trend β of (t) original If the difference between the two exceeds a preset threshold (e.g., |(β)) whitened - β original ) / βoriginal If | > 0.1), it indicates that the pre-whitening process excessively weakens the true long-term trend signal. In this case, a trend-preserving term is introduced into the pre-whitening transformation to correct the segmented pre-whitened sequence. For example, the corrected sequence can be represented as: P wtrend (t)=P w_seg (t)+ α trend ×t, where α trend It is based on β original The trend recovery term is calculated from the autocorrelation coefficient of the sequence (e.g., r1).

[0050] Spatiotemporal registration of meteorological field data from meteorological observations and reanalysis data is performed to obtain a synchronous meteorological element field.

[0051] Specifically, relevant meteorological elements of the study area on one or more isobaric surfaces (e.g., 850 hPa, 700 hPa, 500 hPa) are extracted from reanalysis data (e.g., ERA5). These elements include equivalent potential temperature θe, whole-layer precipitable water volume (IWV), specific humidity q, and zonal and meridional winds u and v. In a preferred embodiment, a dual-track data stream is maintained during spatiotemporal registration to address different data format (point or field) requirements in subsequent calculations: Track 1 (station sequence) uses bilinear interpolation and other methods to convert the reanalysis meteorological elements (e.g., θe, IWV, q) at the grid points. 850 Interpolated to the precise geographical locations of each meteorological observation station, the precipitation sequence P after quality control was obtained. qc (t) A time-synchronized sequence of meteorological elements from the stations, M(t). This sequence M(t) will be primarily used for local thermal and water vapor parameter calculations and change point detection. Orbit 2 (grid data) simultaneously retains the original ERA5 grid data from the surrounding area of ​​the stations (e.g., 5x5 or 3x3 windows). This grid data stream is used to calculate spatial derivatives (such as divergence and vorticity), avoiding the physical and logical contradictions caused by calculating the spatial field from the interpolated station data. External climate oscillation index data are acquired and integrated. For example, monthly values ​​of Nino3.4, PDO, and AMO published by organizations such as NOAA are acquired. When registering these with hourly or daily precipitation sequences, simple linear interpolation may introduce artificial smoothing effects or spurious high-frequency signals. Therefore, in a preferred embodiment, cubic spline interpolation or LOESS smoothing methods can be used to refine the monthly data into daily sequences, generating a climate background index sequence C(t). It should be noted that the C(t) obtained through this processing is only used to characterize the low-frequency evolution trend of the large-scale climate background, providing support for constructing the physical prior probability of abrupt change detection. The synchronous meteorological element field is the collection of various spatiotemporal registration data generated above, including but not limited to the station meteorological element sequence M(t), the climate background index sequence C(t), and the retained original gridded data.

[0052] This embodiment obtains high-quality standardized precipitation time series and synchronous meteorological element fields through quality control, spatiotemporal registration (including dual-track data streams and optimized interpolation methods), and segmented pre-whitening processing with trend preservation verification, providing a reliable data foundation.

[0053] like Figure 2 As shown, according to one aspect of this application, a rainstorm event sequence set is generated, including:

[0054] Based on the thermal and water vapor parameters in the synchronous meteorological element field, a comprehensive meteorological similarity time series is constructed.

[0055] In this embodiment, multiple physical quantities (thermal and water vapor) characterizing atmospheric state are integrated into a single index. Changes in this index can reflect the evolution or passage of weather systems. Specifically, from the station meteorological element sequence M(t), the equivalent potential temperature change Δθe(t), the change in total precipitable water ΔIWV(t), and the change in 850 hPa specific humidity Δq at adjacent times (e.g., t and t-1) are extracted. 850 (t). To enable the synthesis of physical quantities with different dimensions, each sequence of changes is standardized (e.g., divided by its respective standard deviation σ). θe , σ IWV , σ q850 The comprehensive meteorological similarity time series S(t) is constructed by weighted combination. For example, this combination can be a linear function, such as: S(t) = w1 × (Δθe(t) / σ) θe )+w2×(ΔIWV(t) / σ IWV )+w3×(Δq 850 (t) / σ q850 ); where w1, w2, and w3 are weighting coefficients, and their sum is 1, for example, w1=0.3, w2=0.4, and w3=0.3. w2 (total precipitable water content) is given a higher weight because it makes the largest direct contribution to precipitation. Drastic changes in S(t) (i.e., large values ​​or abrupt changes in S(t)) usually indicate a rapid reconfiguration of atmospheric thermodynamic or water vapor conditions, which may correspond to the beginning of a new weather system or the end of an old weather system.

[0056] By analyzing the comprehensive meteorological similarity time series, candidate moments representing weather system transitions are identified, forming a set of weather system transition moments.

[0057] In this embodiment, the moments of drastic change are automatically and objectively identified from the comprehensive meteorological similarity time series S(t). For example, a change point detection algorithm can be used to analyze the comprehensive meteorological similarity time series.

[0058] like Figure 3As shown, in a preferred implementation, forming a weather system transition time set includes:

[0059] An adaptive penalty parameter is dynamically generated based on the local variance sequence and the coefficient of variation of the local variance sequence of the comprehensive meteorological similarity time series.

[0060] Specifically, a sliding window of fixed length (e.g., 24 hours) is used to traverse the comprehensive meteorological similarity time series S(t), calculate the data variance within each window, and construct a local variance sequence σ. local (t). Calculate this σ local The mean (σ) of the (t) sequence local ) and standard deviation std(σ) local ), and obtain its coefficient of variation CV=std(σ local ) / mean(σ local The larger the coefficient of variation (CV), the greater the variation in volatility of the S(t) sequence across different time periods. Based on this CV value, the standard penalty term of the Pruned Exact Linear Time (PELT) algorithm (e.g., β = 2log(N), where N is the total sequence length) is dynamically adjusted to generate the adaptive penalty parameter β. adaptive For example, the following function can be used: β adaptive = (2log(N)) · (1 + k ·CV); where k is an adjustment coefficient, for example, k=0.5. This ensures that when the comprehensive meteorological similarity time series S(t) has large fluctuations (high CV), the penalty parameter is increased accordingly to avoid the algorithm identifying too many false change points in the noise; and vice versa.

[0061] The integrated meteorological similarity time series and adaptive penalty parameters are input into the Pruned Exact Linear Time (PELT) change point detection algorithm to solve for the optimal segmentation path and obtain the candidate change point time series.

[0062] Specifically, the PELT algorithm solves a problem of minimizing a cost function through dynamic programming. Let C(s,t) be the cost of a segment from time s+1 to time t in the comprehensive meteorological similarity time series S(t) (e.g., calculated based on the negative log-likelihood function). For example, the segmented cost function can be: C(s,t) = (ts)×log(σ) st 2 ) +Σ(S(i) -μ st ) 2 / σ st 2 Where s to t is the segmented interval, μ st and σ stLet be the mean and standard deviation of the segment. The goal of the algorithm is to find a set of split points such that the total cost of all segments ∑C(s) i , t i ) and the total penalty term m·β adaptive The sum is minimized, where m is the number of variable points. This is achieved through the recursive relation F(t) = min s<t {F(s) + C(s+1,t) + β} adaptive By combining this with pruning, the optimal splitting path can be obtained in linear time complexity (O(N)), where F(t) is the minimum total cost from the start time to time t, and F(s) is the minimum total cost from the start time to time s. All the splitting points on this path constitute the candidate change point time sequence T. candidate .

[0063] For any candidate change point in the candidate change point time series, examine the significance of the difference between the mean and distribution of the data segments before and after any candidate change point, evaluate the time interval and data stationarity between adjacent candidate change points, eliminate insignificant or false change points caused by noise, and determine the weather system transition time set.

[0064] In this embodiment, the candidate change point time sequence T output by the PELT algorithm is... candidate Post-processing and refinement are performed to improve reliability. Specifically: for the candidate change point time series T candidate Any variable point t in i Extract the S(t) data segments S before and after (e.g., 48 hours each). before and S after Welch's t-test was used to calculate the significance (p-value) of the difference between the means of the two segments. value The Kolmogorov-Smirnov (KS) test was used to verify the difference in the distribution of the two data segments, and D was obtained. ks Statistics and p ks Value, where D ks The KS statistic represents the maximum difference between the cumulative distribution functions of two samples; p ks The p-value of the KS test represents the observed current or more extreme D under the null hypothesis (the two distributions are identical). ks The probability of a value. A comprehensive significance score, such as the Score, can be constructed. sig =-log 10 (p value )×D ks Remove spurious changes, including: removing insignificant changes and retaining the score. sigSignificance points exceeding a certain threshold (e.g., 2.0). As a preferred implementation, to improve adaptability, this significance threshold (2.0) may not be fixed, but rather based on a score. sig The statistical distribution of the sequence itself is dynamically determined; for example, the 90% or 95% quantile (Q) of the sequence can be selected. 90 Or Q 95 This serves as a dynamic threshold. Changes with excessively close intervals are removed, and the time interval Δt between adjacent significant change points is calculated. i = t i+1 - t i If Δt i If the time frame is less than a certain physical minimum (e.g., 6 hours), then the coefficient of variation (CV) of the S(t) segment between these two points is further calculated. segment If CV segment A value very small (e.g., <0.3) indicates that the segment has a gradual change, and this segmentation is likely a false segmentation caused by noise, and should be removed (e.g., merge the interval and remove t). i or t i+1 Remove no-precipitation-response variables and examine the standardized precipitation time series P before and after each variable (e.g., within 12 hours). std If the precipitation change is extremely small (e.g., ΔP < 5 mm / h) and short in duration (e.g., < 3 hours), then the change point may only be meteorological noise rather than a valid weather system transition and should be discarded. The final output is the weather system transition time set T. change .

[0065] Based on the calculation of water vapor flux divergence using synchronous meteorological element fields, the physical continuity before and after each candidate moment in the weather system transition time set is examined. Candidate moments that do not meet the continuity condition are confirmed, and a physical verification transition point set is generated.

[0066] In this embodiment, the weather system transition time set T change Based solely on statistical changes in the station time series S(t), we introduce a spatial physical field (water vapor flux divergence) to verify whether the change truly corresponds to the formation, dissipation, or transformation of a weather system. If the water vapor convergence (the driving force of precipitation) is continuous in both space and time before and after a change point, then the change point should not be considered a dividing point, but rather a continuation of the same precipitation system.

[0067] like Figure 4 As shown, in a preferred implementation, generating the physical verification post-transformation point set includes:

[0068] Based on multi-layer wind field and specific humidity data in the synchronous meteorological element field, a multi-layer water vapor flux divergence field sequence was calculated.

[0069] Specifically, the retained orbital two-grid data is utilized. From this data, the specific humidity (q), zonal wind (u), and meridional wind (v) for several key isobaric surfaces (e.g., 850 hPa, 700 hPa, 500 hPa) are extracted. For each isobaric surface, the zonal component of the water vapor flux (qv) is calculated. x =q×u and meridional component qv y =q×v. The water vapor flux divergence ▽(qV)=Ψ(qv) is calculated using a numerical difference scheme (e.g., central difference). x ) / Ψx+Ψ(qv y ) / Ψy, where Ψ is the partial derivative, and the spatial derivatives Ψ / Ψx and Ψ / Ψy are calculated based on the latitude and longitude spacing Δx and Δy of the grid points. By performing calculations over all time points, a three-layer (or multi-layer) water vapor flux divergence field sequence is finally generated, for example, D 850 (t), D 700 (t), D 500 (t).

[0070] For any candidate moment in the set of weather system transition moments, extract the multi-layer water vapor flux divergence field sequence before and after any candidate moment, and determine the homology from both spatial and temporal dimensions. The homology determination includes: measuring the spatial correlation of the divergence fields before and after, and measuring the temporal evolution consistency of the vertical structure of the divergence fields before and after.

[0071] Specifically, for the weather system transition time set T change Any candidate time t c To measure spatial correlation: extract t c The divergence field D of the lower layers (e.g., 850 hPa) within a short period of time (e.g., 6 hours) 850 (t c-6 :t c ) and D 850 (t c :t c+6 ). Calculate the mean divergence field (or t) for these two time periods. c-1 and t c+1 The spatial correlation coefficient r of the instantaneous field at time (time) spatial Simultaneously, the spatial offset distance d between the extreme centers of water vapor convergence during these two time periods can be identified and calculated. center Optionally, a comprehensive spatial continuity index C can be constructed. spatial For example, C spatial =r spatial ×exp(-d center / L scale ), where L scale The feature scale length (e.g., 200 km). C spatial A higher value indicates greater spatial continuity. For the weather system transition time set T... changeAny candidate time t c To perform a time evolution consistency check: construct a vector that can characterize the vertical structure of the divergence field, such as V. profile (t)=[D 850 (t), D 700 (t), D 500 (t), (D) 700 -D 850 ), (D 500 -D 700 This vector contains the divergence values ​​of each layer and an approximation of the vertical gradient. Calculate t. c The time evolution rate (i.e., time derivative) of the vertical structure vector dV over the first 6 hours before / dt and t c The evolution rate dV in the last 6 hours after / dt. By calculating the cosine of the angle between these two rate vectors, cos(θ) = (dV before dV after ) / (∣dV before ∣×∣dV after The consistency of the evolutionary trend is measured by cos(θ). If cos(θ) is close to 1 (e.g., >0.7), it indicates that the temporal evolution of the vertical structure is consistent over time. c The sequence is continuous, with no mutations occurring.

[0072] If the result of the homology determination is discontinuous, then the candidate moment is confirmed as a valid system transition point and retained in the transition point set after physical verification.

[0073] In this embodiment, for any candidate time t c Based on the overall homology determination results: if the spatial continuity index C spatial The cosine of the included angle, cos(θ), is very high (e.g., >0.7) and t is also very high (e.g., >0.7). c Low-level divergence field (D) before and after (e.g., within 6 hours) 850 All remain negative (i.e., persistent water vapor convergence, for example, <-5×10). -7 kg / (m 2 ·s). When all three conditions above are met simultaneously, then the t is determined to be... c The time before and after belongs to the continuation of the same precipitation system, and is physically continuous. At this time, the t c As a statistical variable point, its physical meaning is negated, and it should be removed from the weather system transition time set T. change Remove from the middle. Conversely, if at least one of the above three conditions is not met (e.g., low spatial field correlation, abrupt change in vertical structure evolution, or interruption of water vapor convergence turning into divergence), then the t is determined to be invalid. cIf the physical continuity of a time point is not met (i.e., it is discontinuous), the candidate time point is identified as a valid system transition point and should be retained. As a preferred implementation, to improve adaptability to different regional climatic characteristics, the fixed threshold (e.g., C) used in the above determination... spatial > 0.7, cos(θ) > 0.7, divergence value < -5×10 -7 All of these can be replaced with dynamic or adaptive thresholds. For example, these thresholds can be obtained based on historical climatological data of the study area, such as using a specific quantile (e.g., the 90th quantile) as the threshold, or introducing a seasonal adjustment factor (e.g., f). season Traverse the weather system transition time set T. change After performing the above checks and screenings on all candidate moments, the set of points that are ultimately retained constitutes the physical verification-based conversion point set T. verified At every moment in this set of points, both the statistical significance test and the physical continuity rejection test are passed.

[0074] Based on the intensity gradient of the standardized precipitation time series and the movement speed of the weather system estimated from the synchronous meteorological field, a dynamic minimum interval duration is constructed.

[0075] In this embodiment, a dynamic time window mechanism is used to make the no-precipitation interval of the segmented rainstorm event no longer a fixed value (such as the traditional 6 hours or 24 hours), but adaptively change according to the current weather system characteristics (precipitation gradient, system movement speed, water vapor conditions).

[0076] In one exemplary embodiment, constructing the dynamic minimum interval duration includes:

[0077] Based on standardized precipitation time series, a multi-scale precipitation gradient sequence reflecting the instantaneous and trend changes in precipitation intensity is generated by using multi-time scale difference combination.

[0078] Specifically, reading the standardized precipitation time series P std (t). Calculate the precipitation intensity gradient at multiple time scales (e.g., 1 hour, 3 hours, 6 hours). For example: 1-hour scale (reflecting instantaneous changes): use forward differencing grad... 1h (t) = P std (t+1) - P std (t); 3-hour timescale (reflecting short-term trends): using central difference grad... 3h (t) = (P std (t+3) - P std (t-3)) / 6; 6-hour timescale (reflecting mid-term trend): using central difference grad 6h(t) = (P std (t+6)- P std (t-6)) / 12; Construct a multi-scale precipitation gradient sequence G by weighted combination precip (t), for example: G precip (t) = 0.5 * |grad 1h (t)| + 0.3 * |grad 3h (t)| + 0.2 * |grad 6h (t)|, where grad 1h (t) represents the precipitation intensity gradient on a 1-hour timescale, grad 3h (t) represents the precipitation intensity gradient on a 3-hour timescale, grad 6h (t) represents the precipitation intensity gradient on a 6-hour timescale. This G... precip The larger the (t) value, the more drastic the change in precipitation intensity.

[0079] Based on multi-layer wind field data in the synchronous meteorological element field, the system velocity sequence is estimated by tracking the movement of the extreme center of water vapor flux divergence.

[0080] In this embodiment, multiple wind fields (such as wind fields u and v at 850 hPa, 700 hPa, and 500 hPa) are extracted from the station meteorological element sequence M(t). An environmental steering airflow can be calculated as a reference, for example: U steer = 0.3*u 850 + 0.5*u 700 + 0.2*u 500 (and V) steer ), where U steer V represents the zonal component of the environmental steering airflow obtained from multi-layer wind field weighted calculations. steer The meridional component of the environmental steering airflow is obtained based on multi-layer wind field weighted calculation. 850 u represents the zonal wind speed component on the 850 hPa isobaric surface. 700 u represents the zonal wind speed component on the 700 hPa isobaric surface. 500 This represents the zonal wind speed component on the 500 hPa isobaric surface. More precisely, this is achieved using the water vapor flux divergence field sequence D. 850 (t), tracking D through pattern matching or correlation methods. 850 The spatial displacement (Δx, Δy) of the extreme center of water vapor convergence (negative divergence center) in (t) between consecutive time points (e.g., time t and time t+3h). The system velocity sequence V is calculated. sys (t), i.e., V sys (t)=sqrt((Δx / 3h) 2 +(Δy / 3h) 2).

[0081] By integrating multi-scale precipitation gradient sequences, system movement velocity sequences, and whole-layer precipitable water extracted from synchronous meteorological element fields, a dynamic minimum interval duration that varies with time is constructed through a function that includes gradient adjustment, movement velocity adjustment, and water vapor adjustment factors.

[0082] Specifically, the basic interval duration T is set. base (e.g., 12 hours). Construct multiple adjustment factors: gradient adjustment factor f grad =exp(-k1*G precip_norm (t)), for example, k1=0.1, where G precip_norm For the standardized precipitation gradient; when the precipitation gradient G precip When the gradient is larger (due to drastic changes in precipitation), the gradient adjustment factor f grad The smaller the value, the smaller the minimum allowable interval T. min The shorter the value, the more sensitively events can be segmented. Speed ​​adjustment factor f speed =1+k2*V sys_norm (t), for example, k2=0.05, where V sys_norm To standardize the system's moving speed; when the system's moving speed V sys The faster the movement speed adjustment factor f is, the higher the speed will be. speed The larger the value, the smaller the minimum allowable interval T. min The longer the interval, the less likely it is to be considered the end of the event, because brief intervals in precipitation caused by the rapid passage of the system should not be regarded as the end of the event. Water vapor moderating factor f vapor =1-k3*IWV norm (t), for example, k3=0.3, where IWV norm This is the standardized whole-layer precipitable water extracted from M(t); when the atmospheric water vapor content IWV is higher (more favorable background conditions), the water vapor moderating factor f vapor The smaller the value, the smaller the minimum allowable interval T. min The shorter the interval, the more likely it is that even a brief pause can be triggered by new convection. It should be noted that, to ensure dimensional consistency in the above formulas, when calculating f... grad f speed f vapor The G used at that time precip (t), V sys Both IWV(t) and IWV(t) preferably use standardized or normalized dimensionless values ​​(e.g., subtracting the mean and dividing by the standard deviation, or dividing by a reference value). In this case, k1, k2, and k3 are dimensionless adjustment coefficients. Finally, a dynamic interval threshold sequence T is constructed. min (t) is: T min (t)=T base ×f grad ×f speed ×fvapor .

[0083] After combining physical verification with the conversion point set, the standardized precipitation time series is segmented to divide the rainstorm event sequence set.

[0084] Specifically, from the standardized precipitation time series P std The traversal begins from the starting point of (t). When precipitation P is detected... std When (t) > 0 (or greater than a certain small threshold), it is marked as the starting point of a rainstorm event. The event is continuously tracked until any of the following conditions are met, at which point it is marked as the ending point: Dynamic interval constraint: A continuous period of no precipitation (or precipitation below a threshold) is detected, and this period exceeds the current dynamic interval threshold sequence T. min (t); Physical change point constraint: The current time t of the traversal reaches the physical verification set T. verified Any time point included in the data; meteorological abrupt change constraint: a drastic change is detected in the comprehensive meteorological similarity time series S(t), for example, |S(t) - S(t-1)| > 2σ. S (where σ) S (where S(t) is the standard deviation). After an event is marked as ended, information such as the start and end times of the event is recorded, and the process continues to iterate from the next moment to identify subsequent rainstorm events. As a preferred implementation, the fixed threshold in the above constraints (e.g., 2σ) S This may lead to statistical consistency issues. Therefore, the threshold can be determined dynamically, for example, by using a peak over threshold (POT) method, or by selecting the 95th percentile (Q) of the |S(t) - S(t-1)| sequence. 95 The adaptive threshold is used as the threshold. Through the above traversal and segmentation, the final rainstorm event sequence set E = {E1, E2, ..., E...} is generated. m}

[0085] Furthermore, attribute features can be extracted from each event in the rainstorm event sequence set E to construct an event attribute feature table F. event .

[0086] Specifically, iterate through each event E in the rainstorm event sequence set E. i From the standardized precipitation time series P std The start and end times t are counted in (t). start and t end Total precipitation R total Maximum hourly rainfall intensity I max Duration D=t end -t startSimultaneously, the event E is extracted from the meteorological element sequence M(t) of the station. i Physical characteristics during the occurrence, such as the average equivalent potential temperature θe mean Maximum total floor precipitation (IWV) max Etc. Integrate all the above attributes to construct the event attribute feature table F. event .

[0087] In an optional embodiment, generating a high-confidence mutation point set specifically includes:

[0088] For a set of rainstorm event sequences, multiple preset mutation detection methods are run in parallel to generate multi-source detection evidence for preliminary mutation candidate points.

[0089] In this embodiment, years that may experience mutations are initially screened from a statistical perspective. Specifically, this is based on the real event attribute feature table F. event Construct one or more annual time series (e.g., annual cumulative heavy rainfall series, or annual heavy rainfall frequency series). For this annual time series, run multiple mature mutation detection methods in parallel (i.e., independently), such as the Mann-Kendall (MK) mutation test, to generate an MK candidate mutation point set C. mk The Pettitt nonparametric test yields the Pettitt candidate point set C. pet Cumulative sum (CUSUM) test to generate CUSUM candidate point set C cusum The Regime Shift test produces a candidate point set C for RS (Regional State Transition). rs The candidate points detected by all the above methods are summarized to form a preliminary complete set of mutation candidate points C. all Preferably, the entire set of preliminary mutation candidate points C... all For each candidate point in the dataset, a mutation intensity index is quantitatively calculated. The magnitude of change for each statistical mutation point is quantified. Specifically, for the entire preliminary set of mutation candidate points C... all Each candidate point t in c Extract precipitation data before and after the mutation (e.g., 5 years each). Calculate the mean μ before the mutation. before and the mean μ after mutation after and the combined standard deviation σ pooled Constructing the mutation intensity index SI(t) c For example, SI(t) c )=∣μ before -μ after ∣ / σ pooled ×(1-p value ), where p value SI(t) represents the minimum significance level calculated for this point using various test methods. cThe larger the SI(t) value, the higher the statistical magnitude and significance of the mutation. All SI(t) values... c The values ​​constitute the mutation intensity vector SI. vector .

[0090] Based on the characteristics of large-scale climate background changes in the climate oscillation index data, we construct the physical prior probabilities of preliminary abrupt change candidate points.

[0091] In this embodiment, the evaluation is performed at t c The likelihood of a real physical mutation occurring at any given moment is entirely independent of statistical evidence.

[0092] In a preferred embodiment, constructing the physical prior probabilities of preliminary mutation candidate points includes:

[0093] For any preliminary abrupt change candidate point, the changes and trend rates of various preset climate indices before and after the preliminary abrupt change candidate point are extracted from the climate oscillation index data.

[0094] Specifically, the climate background index sequence C(t) is read, and three (or more) index sequences, namely Nino3.4, PDO, and AMO, are extracted from it. The initial set of candidate mutation points C is then analyzed. all Any candidate point t in c Calculate the mean exponent ENSO before the mutation (e.g., 24 months). before PDO before AMO before and the mean ENSO after mutation (e.g., 24 months) after PDO after AMO after Therefore, the change ΔENSO = ENSO after -ENSO before (And ΔPDO, ΔAMO). Optionally, t can also be calculated via linear regression. c The slopes of the trends before and after are calculated, and the slope difference Δk is obtained. ENSO Δk PDO Δk AMO All the above changes together constitute the climate index change characteristic matrix M. climate (t c ).

[0095] Based on the change amount and trend change rate, the single index influence intensity of multiple climate indices is quantified, and the coupling effect between the single index influence intensity of multiple climate indices and between multiple indices is calculated to obtain the coupling intensity index.

[0096] In this embodiment, changes in multiple indices are merged into a single driving strength index. Specifically, the influence strength of a single index is calculated, for example, I. ENSO=∣ΔENSO∣ / σ ENSO (where σ) ENSO (This is the historical standard deviation of the index). Calculate pairwise coupling effects, such as I... EP = I ENSO ·I PDO ·sign(ΔENSO·ΔPDO); where the sign function is used to characterize the cooperative effect (when I changes in the same direction). EP (positive) or antagonistic effect (inverse change I) EP (negative). Optionally, calculate the three-exponential coupling I. EPA =(I ENSO ×I PDO ×I AMO (1 / 3)×sign(ΔENSO×ΔPDO×ΔAMO), where ΔENSO is the change in the ENSO index, representing the candidate point t for climate abrupt change. c The difference between the mean of the Nino3.4 index and the mean of the index 24 months prior and 24 months later; ΔPDO is the change in the Pacific Decadal Oscillation Index; ΔAMO is the change in the Atlantic Decadal Oscillation Index; I EP The strength of the ENSO-PDO coupling effect; I ENSO The intensity of the ENSO single-index influence; I PDO The intensity of the single-exponential effect of PDO; I AMO The influence strength of the AMO single index is then determined. Finally, the coupling strength index I is obtained by combining the above factors (e.g., through weighted summation or taking the maximum value). couple (t c ).

[0097] The coupling strength index is mapped to an initial probability value through a piecewise nonlinear function, which is used to characterize the nonlinear influence of coupling effects of different strengths on the mutation probability.

[0098] In this embodiment, the physical driving intensity I couple This is converted to a probability value within the interval [0, 1]. In an exemplary embodiment, the coupling strength index is mapped to an initial probability value using a piecewise nonlinear function, including: adaptively selecting different functions for mapping based on the interval in which the coupling strength index is located. These functions include: a linear function characterizing the influence of weak coupling strength, a sigmoid function characterizing the rapid growth of the influence of moderate coupling strength, and a saturation function characterizing the asymptotic saturation of the influence of strong coupling strength. Specifically, different functions are adaptively selected for mapping based on the interval in which the coupling strength index is located, for example: when I... couple When the value is small (e.g., < 0.5, indicating a weak effect), use a linear segment: P linear =0.2×I couple , where Plinear For the probability value of a linear segment; when I couple Medium (e.g., 0.5≤I) couple When <2.0 (indicating moderate influence), use the sigmoid segment: P sigmoid =0.1+0.6 / (1+exp(-2×(I couple -1.25), this function grows rapidly in the medium influence region, where P sigmoid The probability value of the Sigmoid segment; when I couple When the value is large (e.g., ≥2.0, indicating a strong influence), use the saturation segment: P saturate =0.7+0.1×tanh(I couple -2.0), this function characterizes the asymptotic saturation effect under strong influence, where P saturate This represents the probability value for the saturation segment. Simultaneously, based on whether the change in a pre-defined dominant climate index (e.g., ENSO) exceeds a specific threshold, additional weights are applied to the mapped function value to comprehensively obtain the initial probability value. For example, if |ΔENSO| > 0.5 (indicating significant ENSO activity), an additional weight term P is added. ENSO =0.1×(1-exp(-∣ΔENSO∣)). Combining the above piecewise function values ​​and additional weight terms (e.g., summing them), we obtain the nonlinear prior probability P. nonlinear (t c ).

[0099] Based on the lag correlation analysis between climate oscillation index data and preliminary abrupt change candidate points, the initial probability values ​​are corrected for time lag effects to determine the physical prior probability.

[0100] In this embodiment, considering that the impact of large-scale climate oscillations on regional precipitation often has a lag of several months to several years, specifically, the correlation coefficients between each index in the climate background index sequence C(t) and the regional precipitation sequence at different lag times (e.g., 3, 6, 9, 12 months) are calculated to construct the lag correlation vector R. lag Determine the optimal lag time τ. opt R lag The lag time with the largest absolute value of the correlation coefficient. Using t... c -τ opt Climate index at time (not t) c At time ( ), the initial probability value is recalculated to obtain the time-delay corrected nonlinear prior probability P. nonlinear (t c ), where τ opt This represents the optimal lag time. Optionally, a time delay correction factor λ can also be applied. lag = 0.8 + 0.2·|r(τ opt )|, where r(τ)opt The correlation coefficient is the one corresponding to the optimal lag time. The mutation prior probability vector P is obtained. prior = P nonlinear 修正后 ·λ lag The P prior That is, physical prior information.

[0101] Within the Bayesian posterior probability framework, multi-source detection evidence and physical prior probabilities are integrated to calculate the posterior probability of each preliminary mutation candidate point.

[0102] In a preferred implementation, the posterior probability of each preliminary mutation candidate point is calculated, including:

[0103] Based on the reliability weights and test statistic strengths of different detection methods in multi-source detection evidence, a likelihood function reflecting the support of evidence is dynamically constructed for any preliminary mutation candidate point.

[0104] In this embodiment, the detection results of various statistical methods (MK, Pettitt, etc.) are converted into the likelihood function P(evidence|change) required by the Bayesian formula. Preferably, this includes assigning reliability weights to the detection methods. Specifically, based on historical validation data or expert knowledge, basic reliability weights can be set for four mutation detection methods, for example, the Mann-Kendall method... mK = 0.3, Pettitt method w pet = 0.25, CUSUM method w cusum = 0.25, Regime Shift method w rs = 0.2. Based on each method at candidate point t c The strength of the test statistic at t is dynamically adjusted. For example, if the MK method at t c Point statistics | Z mK If |> 2.58 (i.e., passes the 99% confidence test), then its weight w mK An additional 0.05 can be added; if the Pettitt statistic U pet If the threshold of 95% is exceeded, then w pet Increase by 0.05. After normalization, we obtain t. c Point-specific method weight vector W method (t c Based on this, the likelihood function is dynamically constructed. Specifically, for each candidate mutation point t... c The detection results of the virus in the four methods are statistically analyzed, and a detection vector D(t) consisting of 0-1 indicator variables is constructed. c ) = [d mK d petd cusum d rs ], where d mK To detect indicator variables using the Mann-Kendall method, d pet To detect indicator variables using the Pettitt method, d cusum To detect indicator variables using the CUSUM method, d rs Detect indicator variables using the Regime Shift method. Based on the method weight vector W. method (t c Calculate the weighted detection rate R detect =∑(w i ×d i ), where d i w is the detection indicator variable for the i-th method. i Let be the weight value of the i-th method. Construct the likelihood function P(evidence|change), which can be expressed as the dynamic likelihood value L(t). c For example: L(t) c )=R detect ×exp(-α×(4-∑d i )); where L(t c That is, t c The likelihood function value of a point; α is a penalty parameter (e.g., α=0.3), which reduces the likelihood of points that are only influenced by a minority (e.g., ∑d). i The likelihood value of candidate points detected by the method that is very small; ∑d i Let D(t) be the total number of detection methods. Furthermore, to consider the consistency of the test results, if D(t) c All d in ) i If the mutation direction (enhancing or weakening) indicated by the detection method with a value of 1 is consistent, then the likelihood value L(t) is... c It can be multiplied by an additional reward factor, such as 1.2.

[0105] In Bayesian inference, the likelihood function and physical prior probability are integrated, and the probability distribution is updated through an iterative optimization process until convergence, thus determining the posterior probability of each preliminary mutation candidate point.

[0106] Specifically, read the mutation prior probability vector P prior and dynamic likelihood value L(t) c Applying Bayes' theorem, calculate the posterior probability of the first iteration: P post (1)(t c )=(L(t c )×P prior (t c )) / P(evidence); where P post (1)(t c ) for tc The first posterior probability of a point; P(evidence) is the marginal probability of the evidence, which can be obtained by applying L(t) to all candidate points. c )×P prior (t c The estimation is performed by summing (or weighted averaging). An iterative optimization process is then executed, revising the prior assumptions by utilizing the statistical characteristics of the data itself (e.g., the distribution of time intervals between abrupt changes). Specifically, P... post (1)(t c Points with probabilities greater than a certain intermediate threshold (e.g., 0.5) are considered quasi-mutation points. The temporal distribution characteristics of these quasi-mutation points are analyzed (e.g., the distribution of time intervals between them are calculated), and this distribution characteristic is used to update the mutation prior probability vector P. prior The parameters in the prior distribution are adjusted (e.g., if the actual mutation is found to occur on average once every 10 years, the prior distribution is adjusted to reflect this information). Based on the updated prior probability, Bayes' theorem is applied again to obtain the posterior probability P for the second iteration. post (2). Repeat this iterative process k times until two consecutive iterations P post (k) and P post The set of mutation points generated by (k-1) changes less than the preset convergence threshold (e.g., the change in elements in the set is less than 5%), or the maximum number of iterations is reached (e.g., k=10). The final output is the converged posterior probability distribution P. posterior (t c ).

[0107] Based on the posterior probability distribution, mutation events are screened and confirmed to form a set of high-confidence mutation points.

[0108] As a preferred implementation, an adaptive thresholding method can be used instead of a fixed threshold. Specifically, the natural breakpoint method is used to evaluate the posterior probability distribution P. posterior (t c Cluster analysis is performed on all probability values ​​of the sequence. The split point that minimizes within-group variance and maximizes between-group variance is found, automatically classifying the probability values ​​into three levels: high, medium, and low. From this, a dynamically determined high-confidence threshold θ can be obtained. high (For example, the high and medium level split points calculated by the natural breakpoint method are typically around 0.8) and the medium confidence threshold θ medium (For example, around 0.5). Filter out all P values. posterior (t c )>θ high The points (years) constitute the high-confidence mutation point set T. break In another alternative, simpler implementation, instead of using the natural breakpoint method, a fixed high-confidence threshold can be specified directly, for example, by directly filtering P.posterior (t c Points with a confidence level greater than 0.8 are designated as the high-confidence mutation point set T. break .

[0109] Preferably, after obtaining the high-confidence mutation point set T break Following this, a comprehensive assessment of mutation characteristics is also included. Specifically, this involves integrating the high-confidence mutation point set T. break and mutation strength vector SI vector Construct the mutation intensity assessment matrix M. break This matrix records in detail the comprehensive characteristics of each high-confidence mutation point, including, for example, the mutation time (year), mutation magnitude (from SI). vector The mutation intensity assessment matrix M includes the posterior probability value (confidence level) and the support of the detection method (i.e., which methods detected the point). Further, this mutation intensity assessment matrix M... break It may also include a 95% confidence interval for each mutation point (e.g., estimated using Bootstrap or Monte Carlo methods) and a mutation type label (e.g., based on μ). before and μ after The comparison is labeled as either enhancing or weakening. Mutation strength assessment matrix M break It can be used as input for generating a comprehensive report.

[0110] According to one aspect of this application, a multi-dimensional feature analysis is performed to generate a comprehensive analysis report on rainstorm characteristics, specifically including:

[0111] Perform long-term trend characteristic analysis. This can quantify the overall changing trend of rainstorm characteristics throughout the study period. Specifically, based on the pre-whitening precipitation sequence P... w (t) (this sequence has been de-autocorrelation removed to ensure the validity of the trend test), and a set of rainstorm event sequences E (from which annual rainstorm frequency, annual average intensity, etc. can be extracted). For one or more of the above sequences, the Mann-Kendall (MK) nonparametric test is used to calculate the trend statistic Z-value and significance level p-value. Simultaneously, to quantify the magnitude of the trend, Sen's slope estimation method is used to calculate the rate of change β (e.g., in mm / 10a or times / 10a). A set of trend characteristic parameters T is generated, containing the trend direction (positive or negative Z-value), the magnitude of change (β-value), and the significance level (p-value). trend .

[0112] Perform multi-scale periodic feature decomposition. Identify the oscillation periods at various time scales implicit in the rainstorm sequence. Specifically, for the standardized precipitation time series P stdWavelet analysis is performed on the (t) (or the extracted annual sequence). As a preferred implementation, wavelet packet decomposition (WPD) is employed to obtain a finer frequency division. Exemplary examples include testing various commonly used wavelet bases, such as the Daubechies series (db4, db6), the Symlets series (sym4), etc., before decomposition. The reconstruction error RMSE and signal energy retention rate E are calculated by decomposing and reconstructing the original sequence. ratio Choose (for example) the reconstruction error with the smallest RMSE and E ratio Wavelet bases with a value greater than 0.95 (such as the db4 wavelet base) are used as the optimal wavelet bases for subsequent analysis. optimal Calculate the theoretical maximum number of decomposable layers J. max =floor(log2(N)), where N is the sequence length. Then proceed layer by layer (from j=1 to J). max Perform wavelet packet decomposition and calculate the information entropy H of each layer j. j (Based on the energy percentage p of all nodes in this layer) i Calculate H j =-∑p i ×log(p i When the information entropy gain H j -H j-1 When the value is less than a certain preset threshold (e.g., 0.01), it indicates that the new information that can be provided by continuing the decomposition is limited. In this case, layer j-1 is determined to be the optimal decomposition layer J. opt To reduce boundary distortion at the start and end points of a sequence caused by wavelet decomposition, a symmetric extension method is preferably used at the boundary effects. That is, at P std The two ends of the (t) sequence are each mirrored and extended by N / 10 (or other suitable length) data points. After performing wavelet packet decomposition on the extended sequence, the coefficients corresponding to the extended parts are pruned, and the effective coefficients of the original length are retained to obtain the boundary-corrected wavelet coefficients W. corrected Use the optimal wavelet basis. optimal And the optimal decomposition level J opt J is performed on the (boundary-processed) sequence. opt Layer wavelet packet decomposition yields 2J opt The coefficients of each subband. Calculate the energy E of each subband (j, n). j，n =∑∣W j，n (k)∣ 2 Where n is the sub-band number, W j，n (k) is the subband coefficient sequence. Identify the energy percentage (E). j，n / E total Subbands exceeding a certain threshold (e.g., 5%) are considered significant components, where Etotal The total energy is given. The center period T corresponding to these significant subbands is calculated. j，n The time series of these components can be reconstructed optionally via inverse wavelet packet transform. Finally, the output is a periodic characteristic spectrum S containing the energy distribution of each sub-band. period The dominant periodic sequence P corresponding to the significant components dominant .

[0113] Spatial distribution pattern analysis is performed. The analysis results of individual stations are extended to a surface area to reveal the spatial heterogeneity of rainstorm characteristics. Specifically, the statistical characteristics of each station in the rainstorm event sequence set E (e.g., annual average rainstorm frequency, annual average total precipitation, average intensity, etc.) are interpolated from discrete station locations onto a regular grid (e.g., a 0.1 degree × 0.1 degree grid) of the study area using spatial interpolation methods. Interpolation methods can include, but are not limited to, inverse distance weighted interpolation (IDW), ordinary kriging interpolation, or thin plate spline interpolation. Based on this, the rainstorm frequency, average intensity, and extreme value distribution of each grid point are statistically analyzed. Spatial distribution maps of different levels of rainstorms can be drawn according to preset precipitation thresholds (e.g., distinguishing between different levels of rainstorm events such as 50 mm, 100 mm, and 250 mm based on Chinese meteorological standards). Finally, a spatial distribution feature field G is generated. spatial and intensity level distribution map M intensity .

[0114] Comprehensive feature report generation and output. Specifically, it integrates the trend feature parameter set T. trend Periodic characteristic spectrum S period Spatial distribution characteristic field G spatial and the high-confidence mutation point set T break and mutation intensity assessment matrix M break Compile a time evolution curve (T can be clearly marked on the curve). break (Year of abrupt change), periodic decomposition diagram (e.g., wavelet power spectrum), spatial distribution heatmap (G) spatial Visual charts, etc. Combined with event attribute feature table F event In-depth statistical analysis (e.g., spatiotemporal distribution differences of rainstorm events with different physical causes (e.g., varying θe levels)). Generate a comprehensive analysis report on rainstorm characteristics with rich graphics and clear conclusions. final .

[0115] In a detailed embodiment, it is assumed that a 48-hour standardized precipitation time series P of a certain station is obtained. std (t) (unit: mm / h, standardized so no physical unit, only an example value) and synchronous meteorological element field M(t) (including θe, IWV, q) 850 ) and the corresponding surrounding grid point wind field data. Where Pstd (t) sequence (hourly): t = 0h - 5h: 0; t = 6h - 15h: [0.5, 1.2, 2.5, 3.0, 2.0, 1.5, 1.0, 0.8, 0.4, 0.1] (first precipitation event); t = 16h - 25h: 0 (10-hour interval without precipitation); t = 26h - 35h: [0.3, 1.0, 1.8, 2.2, 1.6, 1.1, 0.7, 0.5, 0.2, 0.1] (second precipitation event); t = 36h - 47h: 0. M(t) sequence: Assume that the comprehensive meteorological similarity time series S(t) calculated from it undergoes drastic changes in statistical characteristics around t=15h (system 1 departs) and around t=25h (system 2 enters). Gridded wind field / specific humidity data (assumptions): t = 9h - 14h: Strong water vapor convergence centers exist (e.g., D). 850 < -7·10 -7 kg / (m 2 •s); t = 16h - 24h: transitions to weak water vapor divergence (e.g., D 850 > 0); t = 26h - 34h: A ​​new water vapor convergence center appears. Assuming the local variance sequence and coefficient of variation (CV) are calculated based on the S(t) sequence, the adaptive penalty parameter β is finally obtained. adaptive = 10.5. Combine the S(t) sequence and β adaptive = 10.5 Input PELT algorithm. Since S(t) undergoes statistical abrupt changes near t=15h and t=25h, the algorithm solves for the candidate change point time sequence T. candidate = {t=15h, t=25h}. Test for t=15h: Extract the S(t) segments before and after this point, and perform Welch's t-test and KS-test. Assume the calculated Score is... sig = 2.8. Meanwhile, P std (t) drops from 0.1 to 0 at this point, indicating a clear precipitation response. This is because 2.8 > 2.0 (or greater than the dynamic threshold Q). 90 The change point is preserved. Test at t=25h: Extract the S(t) segments before and after this point and perform the test. Assume the calculated Score is... sig = 2.5. Meanwhile, P std (t) jumps from 0 to 0.3 at this point, indicating a clear precipitation response. Since 2.5 > 2.0, this change point is retained. The final weather system transition time set T is obtained. change = {t=15h, t=25h}. Perform physical continuity verification to check T. changeDo the statistical changes in the data correspond to actual physical system derailments? Verification at t=15h (the end point of the first precipitation event): Calculate t=9h to t=21h using gridded data (i.e., t... c The water vapor flux divergence field sequence D (± 6h) 850 (t). In the spatial dimension: compare the divergence fields at t=14h (strong convergence field) and t=16h (weak divergence field). Their spatial forms are completely different, and the spatial correlation coefficient r is calculated. spatial Extremely low, for example, r spatial = 0.15; In the time dimension: compare the vertical structure vector V of the two time periods t=9-14h and t=16-21h. profile The evolution rate of (t). The preceding period shows enhanced convergence, while the following period shows divergence; their evolution trends are opposite, resulting in cos(θ) = -0.8 (much less than 0.7). Condition 1 (spatial continuity) is C. spatial (Based on r) spatial =0.15) is much less than 0.7; Condition 2 (time consistency) is cos(θ) = -0.8, which is much less than 0.7; Condition 3 (continuous convergence) is that the water vapor convergence is interrupted after t=15h. The judgment result is discontinuous. The change point t=15h is a valid physical transition point and should be retained. Verification at t=25h (the starting point of the second precipitation event): Calculate from t=19h to t=31h (t c D (± 6h) 850 (t). In the spatial dimension: compare t=24h (weak divergence field) and t=26h (newly generated strong convergence field). Their spatial forms are different, and r is calculated accordingly. spatial = 0.10; In the time dimension: compare the two time periods t=19-24h and t=26-31h. The evolution trend changes abruptly from divergence to convergence, with very low cos(θ), for example, cos(θ) = -0.5. The result is determined to be discontinuous. The change point t=25h is a valid physical transition point and should be retained. Finally, the set of transition points T after physical verification is obtained. verified = {t=15h, t=25h}. Perform dynamic interval determination and event segmentation, and calculate the multi-scale precipitation gradient G. precip (t). This value is positive during the precipitation period (t=6-15, t=26-35) and zero during the intermittent period (t=16-25). Calculate the system's moving speed V. sys (t). Assume a constant velocity of 5 m / s during precipitation periods and 0 m / s during intermittent periods. Extract the total precipitable water volume IWV(t). Assume the standardized value during precipitation periods is IWV. norm =1.5, with an interval of IWV norm = -1.0. Set T base= 12 hours. Calculate T for the interval (t=16-25). min :G precip_norm = 0; V sys_norm = 0; IWV norm = -1.0; f grad = exp(-0.1 · 0) = 1.0; f speed = 1 + 0.05 · 0 = 1.0; f vapor = 1 - 0.3 · (-1.0) = 1.3; T min (Intermission period) = 12 · 1.0 · 1.0 · 1.3 = 15.6 hours. Calculate T for the precipitation period (e.g., t=8h). min :G precip_norm Assume it is 2.0; V sys_norm Assuming a value of 1.0 (the standardized value corresponding to 5 m / s); IWV norm = 1.5; f grad = exp(-0.1 · 2.0) = 0.819; f speed = 1 + 0.05· 1.0 = 1.05; f vapor = 1 - 0.3 · 1.5 = 0.55; T min (Rainfall period) = 12 · 0.819 · 1.05 · 0.55 = 5.67 hours. The dynamic interval threshold sequence T is obtained. min (t), where the threshold remains constant at 15.6 hours during the period t=16-25h. When determining the event boundaries, the traversal begins from t=0h. At t=6h, precipitation P is detected. std (t) > 0, marking the start of event E1: E1 start = 6h. t=15h, P std (t) = 0.1. t = 16h, P std (t) = 0. At this point, the system checks the segmentation condition: dynamic interval constraint: continuous no precipitation duration = 1 hour. T min (t=16h) = 15.6h, the condition (1h > 15.6h) is not satisfied; physical change point constraint: the current time t=16h is not T. verified The point in the equation (point at t=15h). The condition is not met. Weather abrupt change constraint: |S(16) - S(15)| is assumed to be very large (> 2σ). S The conditions are met. The event boundary can also be determined as: t=6h, E1 start = 6h. t=15h, is T verified Point in the middle. Triggering condition (2). Marking event E1 ends: E1end = 15h. Start searching for the next event from t=16h. No precipitation from t=16h to t=24h. t=25h is T verified Point in the middle. Triggering condition (2). t=26h, precipitation P is detected. std (t) > 0, marking the start of event E2: E2 start = 26h. t=35h, P std (t) = 0.1. t = 36h, P std (t) = 0. Start calculating the duration of continuous no precipitation. From t=36h to t=47h, the duration of continuous no precipitation is 12h. During this period, T min (t) is assumed to be 15.6h. Condition (1) (12h > 15.6h) is never satisfied. At t=47h, the data ends. Event E2 marks the end: E2 end = 35h (last precipitation point). Based on T verified The segmentation result for {15h, 25h} is: the rainstorm event sequence set E is: E1: t start = 6h,t end = 15h; E2: t start = 26h, t end = 35h. Event Attribute Characteristic Table F event (Assuming from P) std (t) Calculate the total precipitation and extract example values ​​of θe and IWV from M(t) When the event ID is E1, the start time (t) start The value is 6, and the end time is (t). end The total precipitation was 15, the duration (D) was 10 hours, and the total precipitation (R) was 15. total The value is 12.5 (example value), and the maximum time strength (I) is... max The value is 3.0 (t=9h), the average θe is 340K, and the maximum IWV is 55kg / m. 2 When the event ID is E2, the start time (t) start The value is 26, and the end time is (t). end The value was 35, the duration (D) was 10 hours, and the total precipitation (R) was 35. total The value is 9.3 (example value), and the maximum time strength (I) is 9.3. max The mean value is 2.2 (t=29h), the average θe is 335 K, and the maximum IWV is 48 kg / m. 2 .

[0116] In an optional embodiment, assume T verified = {} (i.e., determining that t=15 and t=25 are the same system, without retaining the change point). t=6h, E1 start= 6h. t=16h, no precipitation, duration = 1h. T min (16h) = 15.6h. 1 < 15.6. No division. ... t = 25h, no precipitation, duration = 10h. T min (25h) = 15.6h. 10 < 15.6. No division. t = 26h, precipitation P std (t) = 0.3 > 0. No precipitation duration reset to zero. Event E1 continues. t = 35h, precipitation P std (t) = 0.1. t = 36h, no precipitation, duration = 1h. T min (36h) = 15.6h. 1 < 15.6. ... t = 47h, no precipitation, duration = 12h. 12 < 15.6. Data ends, event E1 is marked as finished: E1 end = 35h. Output E = {E1}, E1: t start =6h,t end =35h. It can be seen that T verified (Physical variable point constraints) play a crucial and independent role in the segmentation event.

[0117] This embodiment identifies two physically discontinuous precipitation processes (t=6-15h and t=26-35h) as two independent rainstorm events E1 and E2, even if the no-rainfall interval (10 hours) between them is less than the dynamically calculated maximum allowable interval (T). min =15.6 hours). This solves the problem of erroneous merging that may be caused by traditional methods that rely solely on fixed intervals or a single dynamic interval.

[0118] This invention abandons the arbitrary standard of fixed no-precipitation intervals and instead adopts a multi-stage physical constraint mechanism: by constructing a comprehensive meteorological similarity time series and combining it with the PELT algorithm, the timing of possible weather system transitions is statistically identified; physical verification is introduced, and the physical continuity before and after these transition points is verified by calculating the water vapor flux divergence field, thereby stitching together precipitation processes belonging to the same weather system and avoiding erroneous segmentation; a dynamic minimum interval coupled with precipitation gradient and system movement velocity replaces the fixed interval, preventing erroneous merging of precipitation from different systems. This ensures that the final set of rainstorm event sequences is continuous and complete in terms of physical origin, guaranteeing higher physical authenticity of the sequences input for abrupt change detection and reducing the possibility of statistical artifacts from the source. Furthermore, a Bayesian fusion framework is constructed, no longer relying on single statistical test evidence, but introducing physical prior probabilities calculated based on large-scale climate oscillation indices (ENSO, etc.) as an independent judgment criterion. Only when a candidate abrupt change point is supported by multiple sources of statistical evidence and physical mechanisms will it be confirmed as a high-confidence abrupt change point. By filtering out spurious mutation points that have only statistical significance but no physical background support, the reliability of mutation detection is improved.

[0119] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details of the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A storm event identification, spatial typing and driving factor analysis method, characterized in that, The method comprises the following steps: acquiring and processing meteorological observation and reanalysis data to obtain a standardized precipitation time series and a synchronous meteorological element field; based on the standardized precipitation time series and the synchronous meteorological element field, a dynamic time window mechanism based on physical constraints is used to adaptively identify a rainstorm event sequence set; based on the rainstorm event sequence set and external acquired climate oscillation index data, a Bayesian fusion framework is used to integrate statistical test evidence and physical prior information to detect the mutation characteristics of the rainstorm event sequence set, and a high-confidence mutation point set is generated; integrating the rainstorm event sequence set, the high-confidence mutation point set and the standardized precipitation time series, multi-dimensional feature analysis is performed to generate a rainstorm feature comprehensive analysis report; wherein the generation of the rainstorm event sequence set comprises: based on the thermal and water vapor parameters in the synchronous meteorological element field, a comprehensive meteorological similarity time series is constructed; the comprehensive meteorological similarity time series is analyzed to identify candidate time points representing weather system transitions, forming a weather system transition time point set; based on the synchronous meteorological element field, the water vapor flux divergence is calculated to verify the physical continuity before and after each candidate time point in the weather system transition time point set, and the candidate time points that do not meet the continuity condition are confirmed to generate a physically verified transition point set; a dynamic minimum interval duration is constructed according to the intensity gradient of the standardized precipitation time series and the weather system moving speed estimated from the synchronous meteorological element field, and the standardized precipitation time series is segmented based on the physically verified transition point set to divide the rainstorm event sequence set.

2. The method of claim 1, wherein, forming the weather system transition time point set comprises: dynamically generating an adaptive penalty parameter according to the local variance sequence and its coefficient of variation of the comprehensive meteorological similarity time series; inputting the comprehensive meteorological similarity time series and the adaptive penalty parameter into the PELT change point detection algorithm to solve the optimal segmentation path and obtain a candidate change point time sequence; for any candidate change point in the candidate change point time sequence, the significance of the mean and distribution difference of the data segments before and after the candidate change point is verified, and the time interval and data stationarity between adjacent candidate change points are evaluated, and false change points that are not significant or caused by noise are removed to determine the weather system transition time point set.

3. The method of claim 1, wherein, generating the physically verified transition point set comprises: based on the multi-layer wind field and specific humidity data in the synchronous meteorological element field, a multi-layer water vapor flux divergence field sequence is calculated; for any candidate time point in the weather system transition time point set, the multi-layer water vapor flux divergence field sequence before and after the candidate time point is extracted, and homology is determined from the spatial and temporal dimensions, the homology determination including: measuring the spatial correlation of the divergence fields before and after, and measuring the consistency of the vertical structure of the divergence fields before and after; if the result of the homology determination is discontinuous, the candidate time point is confirmed as an effective system transition point and retained in the physically verified transition point set.

4. The method of claim 1, wherein, constructing the dynamic minimum interval duration comprises: based on the standardized precipitation time series, a multi-time scale difference combination is used to generate a multi-scale precipitation gradient sequence reflecting the instantaneous and trend changes of precipitation intensity; based on the multi-layer wind field data in the synchronous meteorological element field, the system moving speed sequence is estimated by tracking the movement of the water vapor flux divergence extreme center; The dynamic minimum interval length varying with time is constructed by a function containing gradient adjustment, moving speed adjustment and water vapor adjustment factors, by integrating the multi-scale precipitation gradient sequence, the system moving speed sequence and the whole-layer precipitable water extracted from the synchronous meteorological element field.

5. The method of claim 1, wherein, The high-confidence mutation point set is generated, including: For the rainstorm event sequence set, a plurality of preset mutation test methods are run in parallel to generate multi-source detection evidence of preliminary mutation candidate points; Based on the characteristics of the climate oscillation index data representing large-scale climate background changes, the physical prior probability of the preliminary mutation candidate points is constructed; In the Bayesian posterior probability framework, the multi-source detection evidence and the physical prior probability are integrated to calculate the posterior probability of each preliminary mutation candidate point; According to the posterior probability, the mutation event is screened and confirmed to form the high-confidence mutation point set.

6. The method of claim 5, wherein, The physical prior probability of the preliminary mutation candidate points is constructed, including: For any preliminary mutation candidate point, the change amount and trend change rate of a plurality of preset climate indices before and after the preliminary mutation candidate point are extracted from the climate oscillation index data; the single-index influence strength of the plurality of climate indices is quantified, and the coupling effect between two or more indices is calculated to obtain a coupling strength index; The coupling strength index is mapped to an initial probability value by a piecewise nonlinear function; Based on the lag correlation analysis between the climate oscillation index data and the preliminary mutation candidate points, the initial probability value is corrected for time lag effect to determine the physical prior probability.

7. The method of claim 6, wherein, The coupling strength index is mapped to an initial probability value by a piecewise nonlinear function, including: According to the interval in which the coupling strength index value is located, different functions are adaptively selected for mapping, including: a linear function for representing the influence of weak coupling strength, a sigmoid function for representing the rapid growth of moderate coupling strength, and a saturation function for representing the gradual saturation of strong coupling strength; and According to whether the change amount of the preset dominant climate index exceeds a certain threshold, an additional weight is applied to the mapped function value to comprehensively obtain the initial probability value.

8. The method of claim 5, wherein, The posterior probability of each preliminary mutation candidate point is calculated, including: Based on the reliability weight and test statistic strength of different detection methods in the multi-source detection evidence, a likelihood function reflecting the evidence support degree is dynamically constructed for any preliminary mutation candidate point; In Bayesian inference, the likelihood function and the physical prior probability are integrated, and the probability distribution is updated through an iterative optimization process until convergence to determine the posterior probability of each preliminary mutation candidate point.

9. The method of claim 1, wherein, The standardized precipitation time series and the synchronous meteorological element field are obtained, including: For the hourly precipitation observation sequence in meteorological observation and reanalysis data, a quality control process based on multi-criteria identification and adaptive interpolation is applied to generate the quality-controlled precipitation sequence; The quality-controlled precipitation sequence is subjected to segmented pre-whitening processing, and after incorporating a trend information retention verification mechanism into the segmented pre-whitening processing, it is standardized to obtain the standardized precipitation time series; The reanalysis meteorological field data in meteorological observation and reanalysis data are spatio-temporally registered to obtain the synchronous meteorological element field.

Images