Method and system for identifying spatial propagation sources of meteorological extreme events

Abstract

The application provides a kind of space propagation source identification method and system based on meteorological extreme event, by time series reconstruction to historical meteorological observation data of research area, generate the meteorological data sequence set containing multiple element meteorological index, based on the extreme event determination strategy of preset, meteorological data sequence set is extracted and handled, obtain meteorological extreme event sequence set;It is calculated for space-time correlation intensity, generate the space-time correlation intensity matrix describing the correlation degree between events in different spatial positions, significant test is carried out to space-time correlation intensity matrix, and the correlation intensity element is screened out, and a significant correlation set is constructed;Based on it, the propagation path is analyzed reversely, and the propagation source position and corresponding propagation influence parameter in the spatial dimension of meteorological extreme event are determined, and the space propagation source identification result is generated.The application can improve the efficiency and positioning accuracy of propagation source identification, and avoid the excessive dependence of traditional model-driven method on initial conditions.

Description

Technical Field

[0001] This invention relates to the field of data processing technology, and in particular to a method and system for identifying spatial propagation sources based on extreme meteorological events. Background Technology

[0002] In the field of meteorological disaster monitoring, identifying the spatial propagation sources of extreme meteorological events is a crucial foundation for formulating disaster prevention and mitigation strategies. Its core lies in determining the location of the propagation source by analyzing the spatiotemporal evolution of meteorological events. Existing technologies are mainly based on forward extrapolation from atmospheric dynamics models or on statistical correlation analysis. However, the data preprocessing stage often suffers from input errors due to the spatiotemporal heterogeneity of multi-source data; the event identification stage often uses single-point threshold judgments, ignoring the spatiotemporal continuity of events; correlation analysis lacks significance testing of correlation strength, making it susceptible to interference from random correlations; and the source location process is highly dependent on initial conditions or empirical parameters, resulting in insufficient accuracy and reliability under complex meteorological conditions, making it difficult to meet the actual requirements of refined disaster early warning for accurate propagation source location. Summary of the Invention

[0003] In view of this, the present invention provides a method and system for identifying spatial propagation sources based on extreme meteorological events. The technical solution of the embodiments of the present invention is implemented as follows:

[0004] On one hand, embodiments of the present invention provide a method for identifying spatial propagation sources based on meteorological extreme events. The method includes: reconstructing historical meteorological observation data of a study area into a time series to generate a meteorological data sequence set containing multiple meteorological indicators, wherein the meteorological data sequence set has a uniform time sampling interval and spatial resolution; performing event extraction processing on the meteorological data sequence set based on a preset extreme event determination strategy to obtain a meteorological extreme event sequence set with spatiotemporal labels, wherein the meteorological extreme event sequence set includes the event occurrence time and corresponding spatial location coordinates; calculating the spatiotemporal correlation strength of the meteorological extreme event sequence set to generate a spatiotemporal correlation strength matrix describing the degree of correlation between events at different spatial locations, wherein the element values ​​of the spatiotemporal correlation strength matrix characterize the event synchronicity level of corresponding spatial location pairs; performing a significance test on the spatiotemporal correlation strength matrix, filtering out correlation strength elements that pass the significance threshold test, and constructing a significant correlation set of meteorological extreme events; and performing reverse propagation path tracing analysis based on the significant correlation set to determine the propagation source location of the meteorological extreme event in the spatial dimension and the corresponding propagation influence parameters, generating a spatial propagation source identification result containing the propagation source coordinates and influence parameters.

[0005] On the other hand, embodiments of the present invention provide a computer system including a memory and a processor, wherein the memory stores a computer program that can run on the processor, and the processor executes the program to implement the steps in the above-described method.

[0006] The spatial propagation source identification method based on meteorological extreme events provided by this invention reconstructs historical meteorological observation data of the study area into a time series, generating a set of meteorological data sequences with a unified time sampling interval and spatial resolution. This eliminates the spatiotemporal heterogeneity caused by inconsistent sampling intervals and spatial resolution differences in multi-source meteorological observation data, providing underlying data consistency assurance for subsequent extreme event extraction and correlation analysis, and improving the reliability of data preprocessing. Based on a preset extreme event determination strategy, the method performs event extraction processing on the meteorological data sequence set, obtaining a set of meteorological extreme event sequences with spatiotemporal labels. Isolated extreme data points are transformed into structured event objects containing the event occurrence time and corresponding spatial location coordinates, establishing a correlation between the event and geographic space. This provides an accurate analysis unit for calculating the spatiotemporal correlation strength, avoiding the analysis errors caused by neglecting spatiotemporal continuity in traditional single-point event identification. The method then performs spatiotemporal correlation strength analysis on the set of meteorological extreme event sequences. The method calculates the spatiotemporal correlation strength matrix, which describes the degree of correlation between events in different spatial locations. It quantifies the correlation between events into matrix element values ​​to characterize the level of event synchronicity, breaking through the limitations of traditional qualitative descriptions or single-dimensional correlation analysis. This enables quantitative modeling of correlation relationships and improves the accuracy of correlation analysis. The spatiotemporal correlation strength matrix is ​​subjected to significance testing, and correlation strength elements that pass the significance threshold test are selected to construct a significant correlation set. This effectively eliminates random correlation interference and ensures that the correlation set has statistically significant reliability, providing a scientific correlation basis for propagation path analysis. Based on the significant correlation set, reverse propagation path tracing analysis is performed to determine the spatial source location of meteorological extreme events and the corresponding propagation influence parameters. This method does not rely on complex physical models and directly infers the propagation source from event correlation data, improving the efficiency and accuracy of propagation source identification and avoiding the excessive reliance on initial conditions in traditional model-driven methods. Attached Figure Description

[0007] Figure 1 This is a schematic diagram illustrating the implementation process of a spatial propagation source identification method based on extreme meteorological events, provided in an embodiment of the present invention.

[0008] Figure 2 This is a schematic diagram of the hardware entity of a computer system provided in an embodiment of the present invention. Detailed Implementation

[0009] This invention provides a method for identifying spatial propagation sources based on extreme meteorological events. This method can be executed by a processor of a computer system. The computer system can refer to devices with data processing capabilities, such as servers, laptops, tablets, and desktop computers.

[0010] Figure 1 This is a schematic diagram illustrating the implementation process of a spatial propagation source identification method based on meteorological extreme events, as provided in an embodiment of the present invention. Figure 1 As shown, the method includes:

[0011] Step S100: Reconstruct the historical meteorological observation data of the study area into a time series to generate a meteorological data sequence set containing multiple meteorological indicators. The meteorological data sequence set has a uniform time sampling interval and spatial resolution.

[0012] Historical meteorological observation data refers to meteorological information collected over a relatively long period within the study area using observational methods such as meteorological stations and meteorological satellites. This data includes observed values ​​of various meteorological elements such as temperature, air pressure, wind speed, and precipitation. Time series reconstruction is the process of reorganizing and processing the original historical meteorological observation data, aiming to ensure that the data has a uniform time sampling interval and spatial resolution. A uniform time sampling interval means that at each fixed point in time, all meteorological indicators have corresponding observational data, ensuring the consistency and comparability of the data in the time dimension, facilitating subsequent time series-based analyses such as trend analysis and periodic analysis. A uniform spatial resolution ensures that meteorological data from different locations have the same measurement standard in terms of spatial scale, enabling comparative and comprehensive analysis of meteorological conditions in different regions within the same framework.

[0013] Step S200: Based on the preset extreme event determination strategy, perform event extraction processing on the meteorological data sequence set to obtain a meteorological extreme event sequence set with spatiotemporal labels. The meteorological extreme event sequence set includes the time of event occurrence and the corresponding spatial location coordinates.

[0014] In one embodiment, step S200 may include the following steps S210 to S260:

[0015] Step S210: Input the meteorological data sequence set into the extreme event detection module, perform threshold crossing detection on each meteorological indicator sequence according to the preset multi-element joint threshold rule, and generate a preliminary extreme event candidate set. The preliminary extreme event candidate set contains meteorological data points that exceed the indicator threshold and their corresponding timestamp information.

[0016] The extreme event detection module is a functional module specifically designed to identify extreme events in meteorological data. It can be a software program or an algorithm model. The preset multi-factor joint threshold rule comprehensively considers the threshold conditions of multiple meteorological indicators. Only when multiple meteorological indicators simultaneously meet their respective threshold requirements is the data point determined to be an extreme event. Threshold crossing detection checks whether each data point in the meteorological indicator sequence exceeds a pre-set threshold. The preliminary extreme event candidate set is selected during the threshold crossing detection process and includes all meteorological data points that exceed the indicator thresholds, along with their corresponding timestamps, which record the time when the extreme event may have occurred.

[0017] In one embodiment, step S210 may include the following steps S211 to S216:

[0018] Step S211: Extract the time series curves of each meteorological indicator from the meteorological data sequence set, and decompose each time series curve using the ensemble empirical mode decomposition method to obtain the intrinsic mode function set and trend component containing different frequency components.

[0019] Time series curves for various meteorological indicators are formed by connecting the observed values ​​of each indicator at different time points, visually demonstrating how the indicator changes over time. Ensemble empirical mode decomposition (EMD) can decompose complex time series curves into multiple intrinsic mode functions (IMFs) with different frequency characteristics and a trend component. IMFs represent the oscillatory components of different frequencies in the time series curve; each IMF has its unique frequency and amplitude, reflecting different fluctuation characteristics in the meteorological data, such as short-term and medium-term fluctuations. The trend component reflects the overall trend of the time series curve, indicating the direction of the meteorological indicator's development over a longer period.

[0020] The time series curves of each meteorological indicator are extracted from the meteorological data series, for example, the time series curve of temperature. Then, the ensemble empirical mode decomposition (EMD) method is used to decompose the curves. Specifically, white noise is first added to the original time series curves to overcome the mode aliasing problem present in traditional EMD methods. After adding white noise, EMD is performed on the series, yielding multiple intrinsic mode functions and one trend component. This process is repeated multiple times, adding different types of white noise each time. Finally, the intrinsic mode functions and trend component obtained from the multiple decompositions are averaged to obtain the final set of intrinsic mode functions and trend component.

[0021] Step S212: Perform sliding window extremum detection on each intrinsic modulus function and trend component, identify the local maxima and local minima of each component within the sliding window, and generate a set of multi-component extremum points.

[0022] Sliding window extremum detection is an effective method for finding local extrema in time series data. It involves setting a fixed-size window and gradually sliding it across the time series, comparing data points within each window to identify local maxima and minima. A multi-component extremum set is a collection of all local maxima and minima detected within the sliding window for each intrinsic modulus function and trend component. For each intrinsic modulus function and trend component, an appropriately sized sliding window is set. Within each window, data points are compared. If the data points within the window first rise and then fall, the last data point in the rising process is the local maximum; if the data points first fall and then rise, the last data point in the falling process is the local minimum. The local maxima and minima detected within all sliding windows for each intrinsic modulus function and trend component are recorded to form the multi-component extremum set.

[0023] Step S213: Compare the extreme points in the multi-component extreme point set with the threshold range of the corresponding component, filter out the abnormal extreme points that exceed the threshold range, and record the timestamp, component type and corresponding meteorological index value of the abnormal extreme points.

[0024] The threshold range for each component is a pre-defined normal value range for each intrinsic modulus function and trend component. Abnormal extreme points refer to those extreme values ​​that exceed the threshold range of their corresponding components. These points often represent anomalies in meteorological data and may be related to extreme weather events. Recording the timestamp, component type, and corresponding meteorological index value of abnormal extreme points is crucial for accurately locating and analyzing these anomalies later.

[0025] Step S214: Perform consistency verification on abnormal extreme points of different components at the same timestamp. When more than a preset number of components have abnormal extreme points at the same timestamp, mark the timestamp as a joint abnormal moment.

[0026] Consistency verification is used to determine whether outlier extreme points of different components at the same timestamp exhibit synergy and correlation. The preset quantity is a pre-defined standard used to measure whether the number of components exhibiting outlier extreme points at the same timestamp reaches a level sufficient to classify it as a joint anomaly. A joint anomaly moment indicates that multiple components simultaneously exhibit outlier extreme points at this time, which may signify a relatively serious anomaly in the meteorological system at that moment, potentially related to the occurrence of extreme weather events. For outlier extreme points of different components at the same timestamp, the number of components exhibiting outlier extreme points is counted. If this number exceeds the preset quantity, the timestamp is marked as a joint anomaly moment.

[0027] Step S215: Perform temporal continuity analysis on the joint abnormal moments, merge consecutive joint abnormal moments into abnormal time periods, and calculate the cumulative deviation value of each meteorological indicator within the abnormal time period.

[0028] Temporal continuity analysis examines whether joint anomalous moments are continuous along the time axis. If multiple joint anomalous moments are temporally adjacent without significant discontinuity, they are merged into a single anomalous period. The cumulative deviation value is the sum of the deviations of the observed values ​​of each meteorological indicator from its normal range within the anomalous period, reflecting the degree to which the meteorological indicator deviates from its normal state during the anomalous period. For joint anomalous moments, they are arranged chronologically to check for continuity between adjacent joint anomalous moments. If continuous, they are integrated into a single anomalous period. For each meteorological indicator within this anomalous period, the deviation of its observed value from its normal range is calculated. For example, if a meteorological indicator has a lower and upper limit for its normal range, for each observed value within the anomalous period, the difference between it and the upper and lower limits of the normal range is calculated, and these differences are then accumulated to obtain the cumulative deviation value of that meteorological indicator within the anomalous period.

[0029] Step S216: Generate an abnormal event score based on the duration and cumulative deviation of the abnormal period, retain abnormal periods whose scores exceed the event score threshold, and generate a preliminary extreme event candidate set containing the time interval and corresponding spatial location of the abnormal period.

[0030] An anomalous event score is a numerical value derived by comprehensively considering the duration and cumulative deviation of the anomalous period, used to assess the severity of the anomalous event. The event score threshold is a pre-set standard used to filter out anomalous periods with higher severity. The preliminary extreme event candidate set consists of anomalous periods whose scores exceed the event score threshold, along with their corresponding spatial location information; these anomalous periods are likely related to extreme weather events.

[0031] Step S220: Perform spatial continuity analysis on the preliminary extreme event candidate set, identify clusters of extreme data points that occur simultaneously in adjacent spatial locations within a preset time window, and generate event units with spatial clustering characteristics. The event unit includes the spatial range and duration interval covered by the cluster.

[0032] Spatial continuity analysis examines whether extreme data points in the initial candidate set of extreme events exhibit a continuous spatial distribution. Adjacent spatial locations refer to geographically close areas. A preset time window is a pre-defined time range used to determine whether extreme data points occur simultaneously within the same time period. An extreme data point cluster is a set of extreme data points occurring simultaneously in adjacent spatial locations within the preset time window. An event unit is composed of extreme data point clusters, exhibiting spatial aggregation characteristics, and includes the specific spatial range covered by the cluster and the duration of the cluster's existence.

[0033] In one embodiment, step S220 may include the following steps S221 to S226:

[0034] Step S221: Construct a spatial grid index structure for the study area, map the anomalous data points in the preliminary extreme event candidate set to the corresponding grid cells, and generate a spatial grid event distribution with time stamps.

[0035] The spatial grid index structure of the study area divides the area into regular grids, each assigned a unique identifier, facilitating rapid spatial location and search. Outlier data points are those exceeding the indicator threshold in the initial extreme event candidate set. Mapping these outlier data points to corresponding grid cells involves assigning them to appropriate grids based on their geographic location information. The time-stamped spatial grid event distribution not only records the presence of outlier data points in each grid cell but also the specific time of their occurrence, allowing for the consideration of both temporal and spatial dimensions in subsequent analysis.

[0036] Step S222: Calculate the global spatial autocorrelation index of the spatial grid event distribution to determine whether there are significant clustering characteristics in the spatial distribution of extreme data points.

[0037] The global spatial autocorrelation index (GSA) is a statistical indicator used to measure whether spatial data distribution exhibits clustering. By analyzing the correlation between grid cells in a spatial grid event distribution, it determines whether extreme data points are spatially random, uniform, or clustered. If the GSA shows significant clustering, it indicates that extreme data points are concentrated in certain areas, potentially suggesting specific meteorological conditions or geographical factors influencing the occurrence of extreme events. Various methods can be used to calculate the GSA, such as Moran's I.

[0038] Step S223: When the global spatial autocorrelation index exceeds the preset clustering threshold, the spatial grid event distribution is clustered using a density-based spatial clustering algorithm. A distance threshold and a density threshold are set, and grid cells with a spatial distance less than the distance threshold and a number of data points exceeding the density threshold are aggregated into an initial extreme data point cluster.

[0039] Density-based spatial clustering algorithms are methods that cluster data points based on their density. A preset clustering threshold is a pre-defined standard used to determine whether the degree of clustering reflected by the global spatial autocorrelation index has reached a level requiring clustering. A distance threshold refers to the maximum allowable distance between two grid cells in space; only when the spatial distance between two grid cells is less than this threshold can they be clustered into the same cluster. A density threshold refers to the minimum number of data points required within a region; only when the number of data points in a grid cell within a region exceeds this threshold will these grid cells be clustered into a single cluster. The initial extreme data point cluster is a preliminary set formed by the clustering algorithm, consisting of grid cells that meet the distance and density threshold conditions.

[0040] When the global spatial autocorrelation index exceeds a preset clustering threshold, a density-based spatial clustering algorithm is initiated. This algorithm traverses each grid cell in the spatial grid event distribution, calculates the spatial distance between adjacent grid cells, and counts the number of data points in each region. If the spatial distance between grid cells in a certain region is less than a distance threshold, and the number of data points in that region exceeds a density threshold, these grid cells are aggregated together to form an initial extreme data point cluster.

[0041] Step S224: Perform boundary identification processing on each initial extreme data point cluster, calculate the spatial contour boundary of the cluster using the convex hull algorithm, extract the grid cell coordinates on the boundary as vertices of the cluster's spatial range, and generate a boundary-optimized extreme data point cluster.

[0042] Boundary identification processing is used to accurately determine the spatial extent of the initial extreme data point cluster. Calculating the spatial contour boundary of the initial extreme data point cluster using the convex hull algorithm yields a concise and accurate representation of the cluster's spatial extent. Extracting the coordinates of the grid cells on the boundary as vertices of the cluster's spatial extent clarifies the specific boundary locations of the cluster; connecting these vertices clearly outlines the cluster's spatial extent. The boundary-optimized extreme data point cluster is obtained after boundary identification processing based on the initial extreme data point cluster, resulting in a more defined and precise spatial extent.

[0043] For each initial cluster of extreme data points, a convex hull algorithm is used. First, the coordinates of the center points of all grid cells within the cluster are used as the input point set. The convex hull algorithm identifies the points among these that form the smallest convex polygon; these points are the boundary points of the cluster's spatial extent. The coordinates of the grid cells corresponding to these boundary points are extracted and used as the vertices of the cluster's spatial extent. By connecting these vertices, the spatial contour boundary of the cluster can be determined, thus generating a boundary-optimized cluster of extreme data points.

[0044] Step S225: Calculate the spatial centroid coordinates and coverage area parameters of the boundary optimization extreme data point cluster. Combine the temporal distribution characteristics of the abnormal data points in the cluster to generate a characteristic extreme data point cluster containing the start timestamp, end timestamp, spatial centroid coordinates, and coverage area parameters.

[0045] The spatial centroid coordinates are the spatial center position of the boundary optimization extreme data point cluster. They are calculated by weighted averaging of the positions of all grid cells within the cluster, reflecting the cluster's spatial concentration trend. The coverage area parameter refers to the size of the spatial area occupied by the boundary optimization extreme data point cluster, which can be obtained by calculating the area enclosed by its spatial outline boundary. The temporal distribution characteristics of outlier data points within the cluster include information such as the start and end times of their occurrence. Characterized extreme data point clusters are formed based on the boundary optimization extreme data point cluster, comprehensively considering the spatial centroid coordinates, coverage area parameter, and temporal distribution characteristics of outlier data points, thus containing more feature information about the cluster. When calculating the spatial centroid coordinates of the boundary optimization extreme data point cluster, the position coordinates of each grid cell within the cluster are multiplied by the number of outlier data points in that grid cell as a weight. Then, the summation of all grid cell coordinates is divided by the total number of outlier data points to obtain the coordinates of the spatial centroid. For the coverage area parameter, geometric calculation methods can be used to calculate its area based on the polygon defined by the spatial outline boundary of the cluster. By combining the temporal distribution of outlier data points within the cluster, the earliest occurrence of an outlier data point is identified as the start timestamp, and the latest occurrence as the end timestamp. The start timestamp, end timestamp, spatial centroid coordinates, and coverage area parameters are then combined to generate a characteristic cluster of extreme data points.

[0046] Step S226: Perform spatial overlap analysis on the characteristic extreme data point clusters. When the spatial overlap area of ​​two characteristic extreme data point clusters exceeds the overlap threshold, they are determined to be spatially correlated clusters. Then, perform a temporal correlation test on the spatially correlated clusters by comparing the overlap interval length between the start and end timestamps of the clusters. When the overlap interval length exceeds the temporal correlation threshold, the two spatially correlated clusters are merged to generate a spatiotemporally correlated extreme data point cluster.

[0047] Spatial overlap analysis aims to determine the degree of spatial overlap between two characteristic extreme data point clusters. The overlap threshold is a pre-defined proportional standard used to determine whether the spatial overlap of two clusters reaches a level that suggests a spatial correlation. If the proportion of the spatial overlap area of ​​two characteristic extreme data point clusters to the area of ​​the smaller unit exceeds the overlap threshold, they are classified as spatially correlated clusters. Temporal correlation testing, based on the identification of spatially correlated clusters, further examines their temporal correlation. The temporal correlation threshold is another pre-defined proportional standard used to determine whether the temporal overlap of two clusters is sufficiently high. When the length of the overlap interval between the start and end timestamps of two spatially correlated clusters exceeds the temporal correlation threshold as a proportion of the shorter cluster's duration, it indicates a strong temporal and spatial correlation between the two clusters. These clusters are then merged into a single spatiotemporally correlated extreme data point cluster, which more comprehensively reflects the spatiotemporal continuity and correlation of extreme events.

[0048] Step S230: Calculate the spatial influence radius parameter based on the spatial range of the event unit, and calculate the event intensity index in combination with the duration interval. The event intensity index is positively correlated with the spatial influence radius parameter and the duration interval length.

[0049] The spatial extent of an event unit refers to the size of the geographical area covered by the event. The spatial radius of influence is an indicator used to measure the spatial extent of an event's influence; it can be calculated by examining the spatial extent of the event unit, such as by calculating the radius of its circumscribed circle. The duration interval refers to the length of time the event lasts from start to finish. The event intensity index is an indicator that measures the severity of an event by comprehensively considering both the spatial radius of influence and the duration interval. It is positively correlated with both the spatial radius of influence and the duration interval, meaning that the larger the spatial radius of influence and the longer the duration, the higher the event intensity index, indicating a greater impact and severity of the event.

[0050] Step S240: Perform temporal continuity verification on the event units, merge adjacent event units that are temporally continuous and spatially overlapping, and generate a set of meteorological extreme event sequences with unique event identifiers.

[0051] Temporal continuity verification checks whether event units are continuous along the timeline. If two event units are sequentially connected in time without a significant time interval and have overlapping areas in space, it indicates that these two event units may be different stages or parts of the same extreme weather event. Merging these temporally continuous and spatially overlapping adjacent event units makes the representation of extreme weather events more complete and accurate. The set of extreme weather event sequences with unique event identifiers is obtained after merging; each event unit is assigned a unique identifier, which facilitates subsequent management and querying of these events.

[0052] For each event unit, they are arranged chronologically, and adjacent event units are checked for temporal continuity. Simultaneously, they are examined for spatial overlap. If both temporal continuity and spatial overlap are met, these two adjacent event units are merged into a new event unit. This process is repeated until all adjacent event units meeting the conditions have been merged. Each merged event unit is assigned a unique event identifier, forming a set of meteorological extreme event sequences with unique event identifiers.

[0053] Step S250: Add spatiotemporal marker information to each event unit in the meteorological extreme event sequence set. The spatiotemporal marker information includes the event start time, end time, spatial center point coordinates, and event intensity index.

[0054] Spatiotemporal labeling information is used to accurately describe the characteristics of extreme meteorological events in time and space. The start and end times of the event clarify the time range of the event, the coordinates of the spatial center point determine the central location of the event in space, and the event intensity index reflects the severity of the event. Adding spatiotemporal labeling information to each event unit makes the set of extreme meteorological event sequences more complete and detailed, facilitating subsequent analysis and research on extreme events, such as analyzing the event's propagation path and impact range.

[0055] For each event unit in the set of meteorological extreme event sequences, the start and end times of the event are extracted from its relevant data. The spatial center coordinates of the event unit are calculated, which can be obtained through geometric calculations of the spatial range covered by the event unit, such as calculating the centroid coordinates of its spatial boundary. Based on the previously calculated event intensity index, this information is combined and added as spatiotemporal labeling information to the corresponding event unit.

[0056] Step S260: Establish an event attribute index table, associate and store spatiotemporal marker information with the unique identifier of the event unit, and the event attribute index table supports fast event retrieval operations based on spatiotemporal range.

[0057] An event attribute index table is a database table structure used to store and manage information related to event units. It links spatiotemporal marker information with the unique identifier of each event unit, allowing for quick retrieval of the corresponding spatiotemporal marker information using the event's unique identifier, and also enabling rapid searching of related event units based on the time and spatial range within the spatiotemporal marker information. The establishment of this index table improves the efficiency of data querying and retrieval, facilitating the rapid location and analysis of extreme meteorological events. When creating the event attribute index table, the table structure is designed to include fields such as the unique identifier of the event unit, the event start time, the end time, the spatial center point coordinates, and the event intensity index. The spatiotemporal marker information and unique identifier of each event unit in the set of extreme meteorological event sequences are stored according to the table structure, establishing the relationships between them.

[0058] Step S300: Calculate the spatiotemporal correlation strength of the meteorological extreme event sequence set to generate a spatiotemporal correlation strength matrix describing the degree of correlation between events at different spatial locations. The element values ​​of the spatiotemporal correlation strength matrix represent the event synchronization level of the corresponding spatial location pair.

[0059] In one embodiment, step S300 may include the following steps S310~S360:

[0060] Step S310: Extract the spatial location coordinates and time interval information of all event units from the meteorological extreme event sequence set, and construct the event spatiotemporal distribution dataset. The event spatiotemporal distribution dataset contains the event occurrence time series and the corresponding event intensity index for each spatial location.

[0061] The spatial coordinates of an event unit clearly define its specific geographical location, while the time interval information records the time range from the start to the end of the event. The event spatiotemporal distribution dataset integrates the spatial coordinates and time interval information of all event units, forming a dataset containing the event occurrence time series and corresponding event intensity index for each spatial location. The spatial coordinates and time interval information of each event unit are extracted from a set of meteorological extreme event sequences. Each spatial location is treated as a record, and the time intervals of all events occurring at that location are arranged chronologically to form an event occurrence time series. Simultaneously, the event intensity index corresponding to each event is recorded. Combining this information constructs the event spatiotemporal distribution dataset.

[0062] Step S320: Perform multi-scale spatial gridding on the event spatiotemporal distribution dataset, dividing the study area into multiple grid cell levels with different resolutions to form a multi-scale spatial analysis framework. Each grid cell level contains the grid cell at the corresponding resolution and the event statistical features within the cell.

[0063] In one embodiment, step S320 may include the following steps S321 to S326:

[0064] Step S321: Extract digital elevation model data of the study area from the geographic information database, and generate a terrain complexity index by calculating the elevation standard deviation and the mean slope. The terrain complexity index is positively correlated with the elevation standard deviation and the mean slope.

[0065] Digital elevation model (DEM) data is a dataset describing the elevation information of a terrain surface and can be obtained from geographic information databases. The standard deviation of elevation is a statistical indicator that measures the dispersion of elevation data within a study area, reflecting the degree of terrain undulation. The mean slope is the average slope at various locations within the study area, reflecting the degree of terrain inclination. The terrain complexity index is an indicator that comprehensively considers both the standard deviation of elevation and the mean slope to measure the complexity of terrain. It is positively correlated with both the standard deviation of elevation and the mean slope; that is, the larger the standard deviation of elevation and the larger the mean slope, the higher the terrain complexity index, indicating a more complex terrain.

[0066] Step S322: Determine the basic grid resolution parameters based on the terrain complexity index. The higher the terrain complexity index, the higher the basic grid resolution parameters.

[0067] The basic grid resolution parameter is the fundamental scale used to divide the study area into grids. The terrain complexity index reflects the complexity of the terrain in the study area. When the terrain complexity index is high, it indicates that the terrain changes are more dramatic, and a finer grid is needed to accurately describe the distribution of terrain and events. Therefore, the basic grid resolution parameter should also be set higher. When the terrain complexity index is low, the terrain is relatively flat, and a finer grid is not needed. The basic grid resolution parameter can be set lower.

[0068] Step S323: Based on the basic grid resolution parameters, generate multiple resolution levels according to the hierarchical progression to form a resolution sequence. The resolution of adjacent levels increases according to a preset ratio, so that the grid units of each level maintain a nested relationship in terms of spatial coverage.

[0069] The resolution sequence consists of multiple layers with different resolutions, each layer increasing in resolution according to a specific pattern. Starting from the base grid resolution parameters, higher resolution layers are generated sequentially according to a preset ratio. The resolution of adjacent layers increases by a preset ratio, ensuring appropriate differences in resolution between different layers. This allows for both a macroscopic overview and detailed microscopic analysis. The grid cells at each layer maintain a nested relationship in spatial coverage, meaning that lower-resolution layers contain multiple higher-resolution layers, enabling analysis and comparison of the same region at different scales.

[0070] Step S324: Perform layer-by-layer meshing on the study area according to the resolution sequence to generate mesh cell levels of different resolutions. Each mesh cell level contains a set of mesh cells and cell boundary coordinates at the corresponding resolution.

[0071] Based on the previously generated resolution sequence, the study area is sequentially meshed. Each resolution level corresponds to a meshing method, generating a set of mesh cells at that resolution and recording the boundary coordinates of each mesh cell. This results in multiple mesh cell levels with different resolutions, each capable of describing the spatial structure of the study area at varying levels of detail. Following the resolution sequence, the study area is meshed sequentially from low to high resolution. For each resolution level, the study area is divided into multiple mesh cells based on the mesh size of that level. The boundary coordinates of each mesh cell are recorded; these coordinates clearly define the specific location and extent of the mesh cell. The set of mesh cells and the cell boundary coordinates for each resolution level are combined to form the corresponding mesh cell level.

[0072] Step S325: Perform spatial encoding processing on each grid cell level, and assign a unique spatial index code to each grid cell through hierarchical spatial indexing. The spatial index code contains the level identifier and the two-dimensional coordinate offset of the cell within the level.

[0073] Spatial encoding is used to facilitate the management and querying of grid cells. A hierarchical spatial index is an index structure used to organize and store grid cell information, enabling rapid location of each grid cell. A unique spatial index code is a unique code assigned to each grid cell, containing a hierarchy identifier and the cell's two-dimensional coordinate offset within that hierarchy. The hierarchy identifier clarifies the resolution hierarchy to which the grid cell belongs, while the two-dimensional coordinate offset determines the grid cell's specific location within that hierarchy.

[0074] When performing spatial encoding for each grid cell level, a hierarchical spatial indexing method is used. A unique level identifier is assigned to each level, and then the offset of each grid cell within that level is calculated based on its two-dimensional coordinate position within that level. The level identifier and the two-dimensional coordinate offset are combined to form a unique spatial index code.

[0075] Step S326: Map the event units in the event spatiotemporal distribution dataset to each grid unit level according to their spatial location, and statistically analyze the distribution characteristics of event occurrence frequency, average event intensity index and event duration in each grid unit to generate a grid event statistical feature table containing spatial index codes.

[0076] Mapping event units in the spatiotemporal distribution dataset to various grid cell levels according to their spatial location means assigning an event to a corresponding grid cell based on its specific location. Statistical analysis of the event frequency, average event intensity index, and event duration distribution characteristics within each grid cell reveals the occurrence and characteristics of extreme events within that cell. A grid event statistical feature table containing spatial index codes associates the spatial index code of each grid cell with its corresponding event statistical features, creating a table that facilitates querying and analysis.

[0077] Each event unit in the event spatiotemporal distribution dataset is mapped to a grid cell hierarchy of different resolutions based on its spatial coordinates. For each grid cell, the number of events occurring within that cell is counted to obtain the event frequency; the average event intensity index of all events within that cell is calculated to obtain the average event intensity index; and the duration of events within that cell is analyzed to obtain the event duration distribution characteristics. These statistical characteristics are then associated with the spatial index codes of the grid cells to generate a grid event statistical characteristic table.

[0078] Step S330: Under each scale grid of the multi-scale spatial analysis framework, calculate the event time difference matrix between any two grid cells. The element values ​​of the event time difference matrix are statistical characteristic parameters of the absolute difference in the occurrence time of events in the corresponding grid cell pair.

[0079] The multi-scale spatial analysis framework includes grid cell levels of different resolutions. The event time difference matrix is ​​calculated at each grid scale, allowing for the analysis of temporal correlations between events at different spatial scales. The event time difference matrix is ​​a two-dimensional matrix, with rows and columns corresponding to different grid cells. Each element in the matrix represents a statistical characteristic parameter, such as the mean or median, of the absolute difference in the occurrence times of events within the corresponding grid cell pair.

[0080] In one embodiment, step S330 may include the following steps S331 to S336:

[0081] Step S331: Extract the event occurrence time series of the target scale grid cell level in the multi-scale spatial analysis framework from the grid event statistical feature table. The event occurrence time series includes the event start timestamp and corresponding event intensity index of all grid cells at that level.

[0082] The target-scale grid cell level is the resolution level that needs to be analyzed within the multi-scale spatial analysis framework. The event occurrence time series is extracted from the grid event statistical feature table, containing the event start timestamps and corresponding event intensity indices for all grid cells in the target-scale grid cell level. This information reflects the occurrence time and intensity of extreme events within this level of grid cells, providing fundamental data for subsequent time difference calculations. Information for the target-scale grid cell level is selected from the grid event statistical feature table. The event start timestamps for each grid cell within this level are extracted and arranged chronologically to form the event occurrence time series. Simultaneously, the event intensity index corresponding to each event start timestamp is recorded.

[0083] Step S332: Perform time series alignment processing on the event occurrence time series of each grid cell, supplement the event state values ​​of missing time points by time interpolation method, and generate an event state sequence with a uniform time interval. The time interval of the event state sequence is consistent with the event sampling frequency of the grid at this scale.

[0084] Time series alignment is used to ensure that the event occurrence time series of each grid cell has a uniform time interval. Temporal interpolation is a technique used to supplement missing event state values ​​at certain time points, estimating the event state at missing time points based on known event time and state information. The event state sequence is obtained after time series alignment, possessing a uniform time interval consistent with the event sampling frequency of the grid at that scale. This ensures that events in different grid cells can be analyzed at the same time scale. For the event occurrence time series of each grid cell, the event sampling frequency (time interval) of that grid scale is first determined. Then, it is checked whether there are missing time points in the event occurrence time series. If missing time points exist, temporal interpolation is used to supplement them, and the supplemented event states are arranged according to the uniform time interval to generate the event state sequence.

[0085] Step S333: Calculate the cross-correlation function of the event state sequences of any two grid cells, identify the time lag corresponding to the maximum value of the cross-correlation function, and use the time lag as the characteristic time difference parameter of the two grid cells. The characteristic time difference parameter characterizes the propagation time feature of the event between the two grid cells.

[0086] The cross-correlation function is used to measure the correlation between two time series. By calculating the correlation between two time series at different time lags, a function reflecting how the correlation changes with time lag is obtained. Identifying the time lag corresponding to the maximum value of the cross-correlation function is equivalent to finding the time difference when the correlation between the two time series is strongest. The characteristic time difference parameter uses this time lag as an indicator to measure the propagation time characteristic of an event between two grid cells, reflecting the time required for an event to propagate from one grid cell to another. For any two grid cell event state sequences, a conventional cross-correlation calculation method can be used: one sequence is shifted on the time axis, multiplied point-by-point with the other sequence, and the results are summed to obtain the correlation values ​​at different time lags. The maximum value in the cross-correlation function is found, and the time lag corresponding to this maximum value is recorded. This time lag is then used as the characteristic time difference parameter between the two grid cells.

[0087] Step S334: Construct an initial time difference matrix based on the characteristic time difference parameters of all grid cell pairs. The rows and columns of the matrix correspond to different grid cells, and the element values ​​are the characteristic time difference parameters of the corresponding grid cell pairs.

[0088] The initial time difference matrix is ​​a two-dimensional matrix constructed based on the characteristic time difference parameters of all grid cell pairs. The rows and columns of the matrix represent different grid cells, and each element in the matrix represents the characteristic time difference parameter of the corresponding grid cell pair. This matrix comprehensively reflects the time difference relationships between all grid cells.

[0089] Step S335: Perform time scale normalization on the initial time difference matrix, converting the matrix element values ​​into proportions relative to the average duration of events under that scale grid, and generate a normalized time difference matrix.

[0090] Time scale normalization is performed to eliminate the influence of differences in event durations across different grid scales on the time difference matrix, making time differences at different scales comparable. This involves converting matrix element values ​​into proportions relative to the average event duration at that grid scale. This is achieved by dividing the characteristic time difference parameter of each grid cell pair by the average event duration at that grid scale, resulting in a proportional value. The normalized time difference matrix is ​​obtained after time scale normalization, with element values ​​representing proportions relative to the average event duration. This allows for comparison and analysis of time differences across different grid cells at the same time scale.

[0091] Calculate the average duration of events at this scale grid. For each element in the initial time difference matrix, divide it by the average duration of the events to obtain a scale value. Refill these scale values ​​into the matrix to generate the normalized time difference matrix.

[0092] Step S336: Calculate the spatial autocorrelation index of the normalized time difference matrix. When the spatial autocorrelation index is lower than the preset correlation threshold, perform spatial smoothing on the matrix and eliminate local noise interference by the neighborhood weighted average method to generate event time difference matrices at each scale.

[0093] The spatial autocorrelation index measures the spatial correlation of elements in a normalized time difference matrix. If the spatial autocorrelation index is below a preset correlation threshold, it indicates the presence of local noise interference in the matrix, resulting in an insufficiently smooth spatial distribution of element values. Spatial smoothing aims to eliminate this local noise interference, making the matrix element values ​​more continuous and stable in space. The neighborhood weighted average method is a spatial smoothing approach that adjusts the value of each element by weighting the neighboring elements. The event time difference matrices at each scale are obtained after spatial smoothing, thus more accurately reflecting the time difference relationships between different grid cells.

[0094] The spatial autocorrelation index of the normalized time difference matrix can be calculated using methods such as the Moran index. The calculated spatial autocorrelation index is then compared to a preset correlation threshold. If it is lower than the threshold, spatial smoothing is required. A neighborhood-weighted averaging method is used: for each element in the matrix, its neighboring elements are selected, and different weights are assigned based on factors such as their distance from the element. A weighted average is then calculated to obtain a new value for the element. After processing all elements in this way, event time difference matrices at various scales are generated.

[0095] Step S340: Based on the event time difference matrix and event intensity index at each scale, calculate the spatiotemporal correlation coefficient at the corresponding scale. The spatiotemporal correlation coefficient is composed of the synchronicity of event occurrence time and the correlation parameters of event intensity by weighted combination.

[0096] The spatiotemporal correlation coefficient is an index calculated by comprehensively considering the temporal and spatial correlation of events. The event time difference matrix at each scale reflects the differences in event occurrence times between different grid cells, while the event intensity index reflects the severity of the event. The spatiotemporal correlation coefficient can be obtained by weightedly combining the synchronicity of event occurrence times and the correlation of event intensity. At each scale, the spatiotemporal correlation coefficient is calculated based on the event time difference matrix and the event intensity index. For the synchronicity of event occurrence times, the element values ​​in the event time difference matrix can be used as a measure; the smaller the time difference, the more synchronized the events. For the correlation of event intensity, the correlation coefficient between the event intensity indices of different grid cells can be calculated. By assigning certain weights to these two parameters and weighting them, the spatiotemporal correlation coefficient is obtained.

[0097] Step S350: Set fusion weights according to the event density of grid cells at each scale, perform weighted fusion of spatiotemporal correlation coefficients at different scales, and generate cross-scale fused spatiotemporal correlation coefficients.

[0098] Event density in grid cells at each scale refers to the frequency of events occurring within that grid cell at each scale. The fusion weights are set based on the event density of grid cells at each scale; scales with higher event density are likely to have a greater impact on the overall correlation and are therefore assigned higher weights, while scales with lower event density have a relatively smaller impact and are assigned lower weights. Weighted fusion of spatiotemporal correlation coefficients at different scales involves multiplying the spatiotemporal correlation coefficient at each scale by its corresponding fusion weight and then summing the results to obtain the cross-scale fused spatiotemporal correlation coefficient. This cross-scale fused spatiotemporal correlation coefficient can integrate information from different scales and more accurately reflect the spatiotemporal correlation of extreme meteorological events.

[0099] Step S360: Standardize the spatiotemporal correlation coefficients of the cross-scale fusion to generate the final spatiotemporal correlation strength matrix. Each element in the matrix corresponds to the degree of correlation of the spatial location pair identified by the row index and column index.

[0100] Standardization is performed to ensure that the spatiotemporal correlation coefficients from cross-scale fusion have a uniform range of values, facilitating comparison and analysis. The final spatiotemporal correlation strength matrix is ​​obtained after standardization. Each element in the matrix corresponds to the degree of correlation between spatial location pairs identified by both row and column indices, typically ranging from 0 to 1. Values ​​closer to 1 indicate a stronger correlation between extreme events at two spatial locations, while values ​​closer to 0 indicate a weaker correlation. Various methods can be used to standardize the spatiotemporal correlation coefficients from cross-scale fusion, such as the min-max standardization method. This involves calculating the minimum and maximum values ​​of the spatiotemporal correlation coefficients from cross-scale fusion, subtracting the minimum value from each element's value, and then dividing by the difference between the maximum and minimum values ​​to obtain a standardized value between 0 and 1. These standardized values ​​are then re-introduced into the matrix to generate the final spatiotemporal correlation strength matrix.

[0101] Step S400: Perform a significance test on the spatiotemporal correlation strength matrix, select the correlation strength elements that pass the significance threshold test, and construct a significant correlation set of meteorological extreme events.

[0102] For example, a random substitution model can be used to test the significance of the spatiotemporal correlation strength matrix. This method tests the significance of statistical results by generating a random substitution dataset with similar characteristics to the original data to simulate possible outcomes under random conditions. Testing the significance of the spatiotemporal correlation strength matrix involves comparing the element values ​​in the matrix with the results generated by the random substitution model to determine whether these element values ​​are statistically significant. The significance threshold is a pre-defined standard used to judge whether an element value is significant. Elements that pass the significance threshold test are selected to form a significant correlation set of meteorological extreme events. The spatial location pairs represented by these elements are more likely to have real correlations rather than being caused by random factors.

[0103] In one embodiment, step S400 may include the following steps S410~S460:

[0104] Step S410: Construct a random alternative dataset with the same temporal and spatial distribution characteristics as the set of meteorological extreme event sequences. The random alternative dataset is generated by randomly rearranging the original event time series, keeping the frequency of event occurrence and spatial distribution density unchanged but shuffling the time order.

[0105] In one embodiment, step S410 may include the following steps S411 to S416:

[0106] Step S411: Extract the time interval information and spatial location coordinates of all event units from the meteorological extreme event sequence set, calculate the time period characteristics of the event occurrence using the periodogram analysis method, and calculate the spatial density distribution characteristics of the event occurrence using the kernel density estimation method.

[0107] The time interval information of an event unit records the time range from the start to the end of the event, while the spatial location coordinates specify the exact geographical location of the event. Periodogram analysis is a technique used to analyze the periodicity of time series. By performing spectral analysis on the event's time series, it identifies the main periodic components, thus revealing the event's periodic characteristics. Kernel density estimation is a non-parametric method used to estimate data distribution density. By processing the spatial location coordinates of event units, it obtains the spatial density distribution characteristics of the event, i.e., the relative density of events occurring at different spatial locations.

[0108] The temporal interval information and spatial coordinates of all event units are extracted from the set of meteorological extreme event sequences. For the temporal interval information, a periodogram analysis method is used. First, the event occurrence time is converted into time series data, and then a Fourier transform is performed on the time series to calculate its power spectral density. The frequency corresponding to the peak in the power spectral density is the main period of the event occurrence. These main periodic information are extracted to obtain the temporal periodic characteristics of the event occurrence. For the spatial coordinates, a kernel density estimation method is used. A kernel function, such as a Gaussian kernel function, is set with the spatial location of each event unit as the center. The density of the surrounding space is estimated based on the distribution of the kernel function. By superimposing the kernel functions of all event units, the spatial density distribution of the event occurrence in the entire study area is obtained, thus yielding the spatial density distribution characteristics of the event occurrence.

[0109] Step S412: Construct an event time point process model based on time period characteristics. The intensity function of the event time point process model is consistent with the time distribution characteristics of the original event, and can generate a random timestamp sequence that conforms to the periodic characteristics.

[0110] An event time-point process model is used to describe the temporal regularity of events. The intensity function determines the probability of an event occurring at different time points. Based on the previously calculated time periodic characteristics, an event time-point process model is constructed such that its intensity function matches the temporal distribution characteristics of the original events; that is, the event occurrence pattern generated by the model has the same periodicity as the event time distribution in the original data. This model can generate random timestamp sequences that conform to periodic characteristics, and these timestamp sequences can simulate the occurrence time of events under random conditions.

[0111] Step S413: Construct an event space sampling model based on the spatial density distribution characteristics, and generate random spatial coordinates that maintain the original spatial density distribution through the rejection sampling method.

[0112] The event space sampling model is used to generate random spatial coordinates that conform to a preset spatial distribution characteristic. This model is constructed based on the previously calculated spatial density distribution characteristics, with the aim of ensuring that the generated random spatial coordinates reflect the spatial distribution pattern of the original events. Rejection sampling is a random sampling technique that can generate random samples conforming to a given probability density function.

[0113] The probability density function of the event spatial sampling model is determined based on the spatial density distribution characteristics. The study area is considered as a two-dimensional plane, and a probability density function is defined on this plane based on the spatial density distribution. The value of this function represents the relative probability of an event occurring at different spatial locations. Random spatial coordinates are generated using a rejection sampling method. First, a rectangular region containing all possible spatial coordinates is defined within the study area, and random points are uniformly generated within this region. For each generated random point, the value of the point under the probability density function is calculated and compared with a randomly generated threshold. If the probability density function value of the point is greater than the threshold, the point is accepted as a valid random spatial coordinate; otherwise, the point is rejected and regenerated. This process is repeated until a sufficient number of random spatial coordinates are generated.

[0114] Step S414: Perform time rearrangement on the event sequence for each spatial location, keeping the number of events unchanged but redistributing the event occurrence time according to the event time point process model to generate a time rearranged event sequence.

[0115] Temporal rearrangement involves reordering the occurrence times of events while maintaining the same number of events at each spatial location. Based on the previously constructed event timeline process model, the occurrence times of events at each spatial location are reassigned so that the new event occurrence times conform to the time cycle characteristics of event occurrence. For the event sequence at each spatial location, the number of events at that location is first counted. Then, a random timestamp sequence with the same number of events at that location is generated using the event timeline process model. These random timestamps are arranged in ascending order and assigned sequentially to the events at that spatial location, thus completing the temporal rearrangement process.

[0116] Step S415: Perform spatial location perturbation processing on the time rearrangement event sequence, adding random perturbation terms that follow a preset distribution to the original spatial location coordinates to generate a spatial perturbation event sequence.

[0117] Spatial location perturbation is used to further simulate random conditions by making small random changes to the spatial locations in a time-rearranged event sequence. A random perturbation term following a preset distribution is added to the original spatial coordinates; this distribution can be Gaussian, uniform, or similar. The resulting spatially perturbed event sequence exhibits some randomness in spatial location, but overall retains spatial distribution characteristics similar to the original data.

[0118] For each event in the time-rearranged event sequence, its spatial location is represented by the original spatial coordinates. A pre-defined distribution, such as a Gaussian distribution, is chosen, and its mean and standard deviation are determined. For each event's spatial coordinates, a perturbation term randomly drawn from the pre-defined distribution is added to both its x and y axes. The spatial coordinates with the added perturbation term are then used as the new spatial location, thus generating a spatially perturbed event sequence.

[0119] Step S416: Perform spatiotemporal consistency verification on the spatial perturbation event sequence. Use autocorrelation analysis to check whether the autocorrelation characteristics of the time series are consistent with the original sequence. Use the spatial autocorrelation index to check whether the spatial distribution characteristics are maintained. Delete event units that do not meet the consistency requirements and generate the final random substitution dataset.

[0120] Spatiotemporal consistency verification ensures that the generated spatially perturbed event sequence shares similar characteristics with the original data in both time and space. Autocorrelation analysis is used to examine the autocorrelation characteristics of the time series, i.e., the correlation between different time points in the time series. By comparing the time series autocorrelation characteristics of the spatially perturbed event sequence with those of the original sequence, it is determined whether they are consistent in the temporal dimension. The spatial autocorrelation index measures the clustering of spatial data distribution; by comparing the spatial autocorrelation index of the spatially perturbed event sequence with that of the original data, it is determined whether they are consistent in the spatial dimension. Deleting event units that do not meet consistency requirements ensures that the final randomized substitute dataset accurately simulates random conditions and shares similar spatiotemporal characteristics with the original data.

[0121] Step S420: Repeat the spatiotemporal association strength calculation process based on the random substitution dataset to generate multiple random association strength matrices. The number of random association strength matrices is determined by the preset number of repeated samplings.

[0122] Based on the previously generated final random substitution dataset, the same spatiotemporal association strength calculation process as with the original data is repeated. This includes multi-scale spatial gridding of the random substitution dataset, calculating the event time difference matrix, and spatiotemporal association coefficients, ultimately generating a random association strength matrix. The preset number of resampling iterations determines the number of random association strength matrices generated. Through multiple repetitions, the distribution of spatiotemporal association strength under random conditions can be obtained. The final random substitution dataset is processed multiple times according to the preset number of resampling iterations. Each processing step strictly follows the spatiotemporal association strength calculation steps for the original data. First, multi-scale spatial gridding is performed, mapping event units in the random substitution dataset to grid unit levels of different resolutions, and statistically analyzing the event characteristics within each grid unit. Then, the event time difference matrix is ​​calculated, obtaining the event time difference matrices at each scale through time series alignment, cross-correlation function calculation, and other steps. Next, the spatiotemporal association coefficients are calculated, weighted fusion, and standardized to finally generate the random association strength matrix.

[0123] Step S430: Extract the association strength value corresponding to the position of each element in the original spatiotemporal association strength matrix from each random association strength matrix, and construct a random distribution model for each element. The random distribution model includes the probability density function and cumulative distribution function of the association strength value.

[0124] From each random association strength matrix, the association strength value at the corresponding position is extracted according to the position index of the element in the original spatiotemporal association strength matrix. For each element in the original spatiotemporal association strength matrix, the association strength values ​​at the corresponding positions in all random association strength matrices are collected, and these values ​​constitute a random sample set for that element. Based on this random sample set, a random distribution model for each element is constructed. The random distribution model includes a probability density function and a cumulative distribution function. The probability density function describes the probability distribution of the association strength value over different values, and the cumulative distribution function represents the probability that the association strength value is not greater than a certain set value.

[0125] For each element in the original spatiotemporal correlation strength matrix, iterate through all random correlation strength matrices and extract the correlation strength value at the corresponding position. Collect these values ​​to form a random sample set for that element. Based on the characteristics of the random sample set, select an appropriate method to construct a random distribution model. If the sample set conforms to a normal distribution, a parametric method can be used to determine the parameters of the normal distribution by calculating the mean and standard deviation of the samples, thereby obtaining the probability density function and the cumulative distribution function. If the sample set does not conform to a normal distribution, a nonparametric kernel density estimation method is used. A kernel function is set with each sample point as the center. By superimposing and normalizing the kernel functions, the probability density function is obtained, and then the cumulative distribution function is calculated.

[0126] Step S440: Calculate the quantile of each element value in the original spatiotemporal correlation strength matrix in the corresponding random distribution model, and generate a significance p-value. The significance p-value represents the probability that the element value is generated by a random process.

[0127] The p-value is an important indicator used to determine whether the element values ​​in the original spatiotemporal correlation strength matrix are statistically significant. It involves calculating the quantile of each element value in the corresponding random distribution model, thus determining the relative position of that element value within the random distribution. The p-value, derived from the quantiles, represents the probability that the element value was generated by a random process. If the p-value is small, it indicates that the element value is unlikely to be caused by random factors and has high statistical significance.

[0128] In one embodiment, step S440 may include the following steps S441 to S446:

[0129] Step S441: Extract all association strength values ​​corresponding to the same grid cell pair from multiple random association strength matrices, and construct a random sample set of association strengths for the grid cell pair.

[0130] For each grid cell pair in the original spatiotemporal correlation strength matrix, the correlation strength value at the corresponding position is extracted from multiple random correlation strength matrices. These values ​​are collected to form a random sample set of the correlation strength of that grid cell pair. This sample set reflects the possible range of correlation strength values ​​for that grid cell pair under random conditions, providing basic data for subsequent significance testing. The process involves traversing multiple random correlation strength matrices, and for each grid cell pair in the original spatiotemporal correlation strength matrix, extracting the correlation strength value at the corresponding position in all random correlation strength matrices based on its position index.

[0131] Step S442: Perform a normality test on the random sample set of association strength. If the sample set conforms to a normal distribution, construct a normal distribution model using parametric methods; otherwise, construct a probability distribution model using nonparametric kernel density estimation methods.

[0132] Normality testing is used to determine whether a random sample set of association strength conforms to a normal distribution. If it does, parametric methods can be used to determine the parameters of the normal distribution by calculating the sample mean and standard deviation, thus constructing a normal distribution model. If it does not conform to a normal distribution, nonparametric kernel density estimation methods are needed. A kernel function is set centered on each sample point, and the kernel functions are superimposed and normalized to construct a probability distribution model. Normality testing methods, such as the Shapiro-Wilk test, are used to test the random sample set of association strength. If the test results indicate that the sample set conforms to a normal distribution, the sample mean and standard deviation are calculated. A normal distribution model is constructed with the mean as the center and the standard deviation as the scaling parameter. The probability density function of this model can be expressed using the formula for a normal distribution. If the sample set does not conform to a normal distribution, a suitable kernel function, such as the Gaussian kernel function, is selected. Density estimation is performed on the surrounding space centered on each sample point based on the distribution of the kernel function. The kernel functions of all sample points are superimposed and normalized to obtain the probability density function of the probability distribution model.

[0133] Step S443: Input the correlation strength value of the corresponding grid cell pair in the original spatiotemporal correlation strength matrix into the probability distribution model, calculate the probability integral of the right tail of the value, and generate the initial significance p value.

[0134] Substitute the association strength values ​​of the corresponding grid cell pairs in the original spatiotemporal association strength matrix into the previously constructed probability distribution model. Calculate the probability integral of the right-hand side of this value, which is the probability that it is greater than the association strength value. This probability is the initial significance p-value, representing the probability of obtaining a value greater than or equal to the association strength value under random conditions. For each grid cell pair in the original spatiotemporal association strength matrix, input its association strength value into the corresponding probability distribution model. If it is a normal distribution model, calculate the probability of the right-hand side of the value based on the cumulative distribution function of the normal distribution. If it is a probability distribution model constructed using nonparametric kernel density estimation, calculate the probability by integrating the probability density function to the right of the value.

[0135] Step S444: Perform multiple test correction on the initial significance p-values ​​of all grid cell pairs, adjust the significance threshold by controlling the false discovery rate, and generate the corrected significance threshold.

[0136] Because significance testing requires examining the association strength values ​​of multiple grid cell pairs, there is a problem with multiple testing. Multiple testing can lead to an increase in false positives, i.e., incorrectly determining that some association strength values ​​are significant. To control the false discovery rate, the initial significance p-values ​​of all grid cell pairs need to be corrected using multiple testing. False discovery rate control methods, such as the Benjamini-Hochberg method, adjust the significance threshold based on the initial significance p-values ​​and the total number of tests to generate a corrected significance threshold.

[0137] Step S445: Compare the initial significance p-value with the corrected significance threshold. When the initial significance p-value is less than the corrected significance threshold, the association strength element is determined to pass the significance test.

[0138] The initial significance p-value of each grid cell pair is compared with the corrected significance threshold. If the initial significance p-value is less than the corrected significance threshold, it indicates that the association strength element is unlikely to be caused by random factors and has high statistical significance. Therefore, the association strength element is deemed to have passed the significance test.

[0139] Step S446: Record the position index, original association strength value, and corrected significance p-value of the association strength elements that pass the significance test, and generate a list of significance test results.

[0140] The relevant information of the association strength elements that passed the significance test is recorded, including their position index in the original spatiotemporal association strength matrix, the original association strength value, and the corrected significance p-value. This information is then compiled into a list, namely the significance test results list, for convenient subsequent analysis and use.

[0141] For each association strength element that passes the significance test, record its row and column indices in the original spatiotemporal association strength matrix. These two indices together determine the element's position. Also record the element's original association strength value and the corrected significance p-value. Organize this information into a list, such as a table, where each row represents an association strength element that passed the significance test, including its position index, original association strength value, and corrected significance p-value.

[0142] Step S450: Filter out the association strength elements with a significance p-value less than the preset significance threshold, and record the grid cell pairs and association strength values ​​corresponding to these elements.

[0143] The preset significance threshold is a pre-defined standard used to further filter out elements with high association strength. The association strength elements in the significance test result list are compared with the preset significance threshold according to their corrected significance p-values, and elements with significance p-values ​​less than the preset significance threshold are filtered out. The corresponding grid cell pairs and their association strength values ​​are recorded; the associations between the grid cell pairs represented by these elements are more likely to be real and have high significance.

[0144] Step S460: Perform spatial connectivity analysis on the selected association strength elements, construct an association network using graph theory, retain the set of association strength elements that form connected components, delete isolated association strength elements, and construct the final significant association set of meteorological extreme events.

[0145] Spatial connectivity analysis aims to determine whether the selected correlation strength elements are spatially connected. Using graph theory, each grid cell is treated as a node in a graph, and the connections between grid cell pairs corresponding to correlation strength elements are considered as edges, thus constructing a correlation network. A connected component is the largest subgraph composed of interconnected nodes. Retaining the set of correlation strength elements forming a connected component means keeping only those spatially related relationships that form a whole. Isolated correlation strength elements are removed because these elements may be due to random factors and lack actual correlation significance. The final set of significant correlations for extreme meteorological events is obtained after spatial connectivity analysis, containing correlation strength elements with high spatial connectivity and significance.

[0146] The selected association strength elements are converted into nodes and edges in a graph. Each grid cell is treated as a node, and the associations between grid cell pairs represented by the association strength elements are treated as edges. The weight of each edge can be set to the association strength value. Graph theory algorithms, such as Depth-First Search (DFS) or Breadth-First Search (BFS), are used to find the connected components in the graph. The set of association strength elements that form connected components is retained, while isolated association strength elements that do not belong to any connected component are deleted.

[0147] Step S500: Based on the significant association set, perform reverse tracing analysis of the propagation path to determine the spatial source location of the meteorological extreme event and the corresponding propagation influence parameters, and generate spatial propagation source identification results containing the source coordinates and influence parameters.

[0148] In one embodiment, step S500 may include the following steps S510~S570:

[0149] Step S510: Convert the grid cell pairs in the significant association set into a directed graph structure, where nodes represent grid cells, directed edges represent significant associations between grid cell pairs, and edge weights are association strength values.

[0150] A directed graph structure is a graph model used to represent directional relationships between nodes. This method converts grid cell pairs in a salient set of connections into a directed graph structure, where each grid cell corresponds to a node in the graph. The salient relationships between grid cell pairs are represented by directed edges, with the direction of the edge indicating the direction of the connection, and the weight of the edge set being the connection strength value. This allows for a visually intuitive representation of the salient set of connections, facilitating subsequent path analysis and source identification. For each grid cell pair in the salient set of connections, the two corresponding grid cells are treated as two nodes in the directed graph. The direction of the directed edge is determined based on the directionality of the connection, for example, from one grid cell to another. The connection strength value of the grid cell pair is used as the weight of the directed edge.

[0151] Step S520: Based on the time order of events in the meteorological extreme event sequence set, add directional attributes to the edges of the directed graph, pointing from the grid cell with the earlier event time to the grid cell with the later event time, and construct a spatiotemporal propagation directed graph.

[0152] Spatiotemporal propagation directed graphs are constructed by incorporating temporal information about extreme meteorological events into directed graphs. Based on the temporal order of events within each grid cell of the extreme meteorological event sequence set, directional attributes are added to the edges of the directed graph, causing the edges to point from grid cells with earlier event times to those with later event times. This method of constructing a spatiotemporal propagation directed graph more accurately reflects the possible spatial and temporal propagation paths of extreme meteorological events.

[0153] In one embodiment, step S520 may include the following steps S521 to S527:

[0154] Step S521: Extract the event time series of each grid cell from the set of meteorological extreme event sequences, calculate the time difference between the start time of each event cell and the start time of the study period, and generate the event relative time parameters.

[0155] The event time series records the timing information of extreme meteorological events within each grid cell. The time difference between the start time of each event cell and the start time of the study period is calculated to obtain the relative time parameter of the event. This parameter unifies the timing of events onto a relative time scale, facilitating comparisons of the chronological order of events within different grid cells.

[0156] Extract the event time series for each grid cell from the set of meteorological extreme event sequences. These sequences contain the start time of each event within that grid cell. Determine the start time of the study period, for example, using the specific date and time the study begins. For each event cell, calculate the time difference between its start time and the start time of the study period. This time difference can be converted into a value in hours, days, etc., to obtain the relative time parameters of the event.

[0157] Step S522: Perform statistical analysis on the relative time parameters of events for each grid cell, determine the characteristic event time of the grid cell using the median method, and generate a node time attribute vector.

[0158] The median method is a statistical analysis method used to find the median value of a dataset. It involves statistically analyzing the relative time parameters of events for each grid cell and using the median method to determine the characteristic event times for that grid cell. These characteristic event times represent the typical times when events occur within that grid cell. The characteristic event times of each grid cell are combined to generate a node time attribute vector, which provides a basis for subsequently determining the direction of edges in a directed graph. For each grid cell, the relative time parameters of all its event cells are collected to form a dataset. Using the median method, this dataset is arranged in ascending order, and the middle value is identified as the characteristic event time for that grid cell. If the number of data points in the dataset is even, the average of the two middle values ​​is taken as the median. The characteristic event times of each grid cell are then arranged in the order of the grid cells to form a node time attribute vector.

[0159] Step S523: For each pair of grid cells in the significant association set, compare the node time attribute values ​​of the two grid cells. When the node time attribute value of the first grid cell is less than that of the second grid cell, determine the direction as from the first grid cell to the second grid cell.

[0160] For each pair of grid cells in the significantly associated set, compare the node time attribute values ​​of the two grid cells based on the previously generated node time attribute vector. If the node time attribute value of the first grid cell is less than that of the second grid cell, it means that the event in the first grid cell occurred relatively earlier. According to the event propagation logic, determine the direction of the directed edge from the first grid cell to the second grid cell.

[0161] Step S524: Calculate the time delay parameter of the directed edge. The time delay parameter is the ratio of the difference between the time attribute values ​​of two grid cell nodes to the average duration of the event.

[0162] The time delay parameter measures the relative time required for an extreme meteorological event to propagate from one grid cell to another. The time delay parameter for a directed edge is calculated by dividing the difference in the node time attribute values ​​of the two grid cells by the average duration of the event. The average duration of the event is the average of the durations of all extreme meteorological events. By calculating the time delay parameter, the temporal characteristics of event propagation between different grid cells can be described more accurately.

[0163] Step S525: Construct an edge weight calculation model. Input the association strength value and time delay parameter into the model to generate edge weight parameters. The calculation formula for the edge weight parameters is the association strength value multiplied by one and the difference between the time delay parameter and the time delay parameter.

[0164] The edge weight calculation model is used to determine the weight of directed edges by comprehensively considering the association strength value and the time delay parameter. The association strength value represents the tightness of the association between two grid cells, while the time delay parameter represents the relative time it takes for an event to propagate between the two grid cells. The edge weight parameter is calculated by multiplying the association strength value by the difference between the association strength value and the time delay parameter. This setting aims to consider both the association strength and the event propagation time factor. If the time delay parameter is large, it indicates that the event propagation time is long, and the edge weight will be correspondingly reduced; if the time delay parameter is small, it indicates that the event propagation time is short, and the edge weight will be relatively large.

[0165] Step S526: Construct a weighted adjacency matrix for the spatiotemporal propagation directed graph using edge weight parameters. The rows and columns of the weighted adjacency matrix correspond to grid cells, and the element values ​​are edge weight parameters. Grid cell pairs that do not have significant associations have an element value of 0.

[0166] A weighted adjacency matrix is ​​used to represent the connectivity relationships and edge weights between nodes in a graph. Using the edge weight parameters calculated earlier, a weighted adjacency matrix for a spatiotemporally propagating directed graph is constructed. The rows and columns of the matrix correspond to different grid cells, and the element values ​​are the edge weight parameters between corresponding grid cell pairs. If there is no significant association between two grid cells, i.e., no corresponding directed edge, then the element at that position in the matrix has a value of 0.

[0167] Create a square matrix with the same number of grid cells as the weighted adjacency matrix. For each directed edge in the directed graph, fill the corresponding position in the weighted adjacency matrix with the edge weight parameter based on the row and column indices of its corresponding grid cell pair (the row index corresponds to the starting grid cell, and the column index corresponds to the ending grid cell). For grid cell pairs that are not significantly related, set the element value at the corresponding position in the matrix to 0.

[0168] Step S527: Sparsify the weighted adjacency matrix by retaining the top K largest weight edges of each node and deleting edges with smaller weights to control the complexity of the directed graph. K is an adaptive parameter determined based on the average connectivity of the nodes.

[0169] Sparsity reduction aims to reduce the number of non-zero elements in the weighted adjacency matrix, thereby lowering the complexity of the directed graph. Retaining the top K weighted edges of each node means that for each grid cell (node), only the K edges with the highest weights associated with it are kept, while edges with lower weights are deleted. K is an adaptive parameter determined based on the average node connectivity. Average node connectivity refers to the average number of other nodes connected to each node. Sparsity reduction makes the directed graph more concise, highlighting important relationships and facilitating subsequent path discovery and analysis.

[0170] Calculate the average node connectivity, which is the average number of edges connecting all nodes. Determine the adaptive parameter K based on the average node connectivity; for example, K can be set to a multiple of the average node connectivity. For each row (corresponding to one node) in the weighted adjacency matrix, sort the elements of that row in descending order, retain the K largest elements, and set the values ​​of the remaining elements to 0.

[0171] Step S530: Perform path mining on the spatiotemporal propagation directed graph, set a path weight threshold, extract all propagation paths whose cumulative weight exceeds the threshold, and generate a propagation path set.

[0172] Path mining involves identifying potential propagation paths in a directed graph of spatiotemporal propagation. Setting a path weight threshold helps filter out propagation paths of significant importance. The cumulative weight of a propagation path is the sum of the weights of all edges along that path. All propagation paths with cumulative weights exceeding the threshold are extracted and grouped together to form a propagation path set. These propagation paths represent the possible propagation trajectories of extreme meteorological events and are crucial for analyzing the source and process of event propagation.

[0173] In a spatiotemporally propagating directed graph, a graph search algorithm, such as depth-first search or breadth-first search, is used to traverse all possible paths in the graph. For each path, its cumulative weight is calculated, which is the sum of the edge weights of all edges on the path. The cumulative weight is compared with a preset path weight threshold; if the cumulative weight exceeds the threshold, the path is recorded. All propagation paths that meet the criteria are compiled into a set, i.e., the propagation path set.

[0174] Step S540: Perform cluster analysis on the paths in the propagation path set. Use spectral clustering to group propagation paths with similar spatial orientations into one class to generate propagation path clusters. Calculate the average path of each path cluster as the main propagation path.

[0175] Cluster analysis is the process of grouping similar data objects. For paths in a propagation path set, spectral clustering is used to group paths with similar spatial orientations into one class. Spectral clustering is a method of clustering based on the spectral features of a graph. By calculating the similarity matrix between paths, it maps the paths to a low-dimensional space for clustering. The generated propagation path clusters are sets composed of propagation paths with similar spatial orientations. The average path of each path cluster is calculated as the dominant propagation path. The dominant propagation path can represent the typical characteristics of the paths in the cluster, helping to more clearly understand the main propagation direction of extreme meteorological events.

[0176] Each path in the propagation path set is treated as a data object, and a similarity matrix is ​​constructed between the paths. The similarity between paths can be calculated by comparing factors such as the sequence of grid cells traversed by the path and the spatial orientation of the path. Using spectral clustering, the similarity matrix is ​​decomposed into eigenvectors, and the paths are clustered based on these eigenvectors to generate propagation path clusters. For each propagation path cluster, the average path value of all paths in the cluster is calculated. This can be achieved by averaging the coordinates of each node on the path to obtain the position of the node on the average path. Connecting these nodes sequentially yields the main propagation path.

[0177] Step S550: On each main propagation path, trace back from the end point to the starting point using the reverse tracing method, calculate the propagation contribution of each node on the path, taking into account the node's position, edge weight, and number of paths; sort the nodes according to their propagation contribution, and identify the node with the largest propagation contribution as a candidate for the propagation source location.

[0178] The reverse tracing method starts from the end of the main propagation path and traces back to the starting point, analyzing the contribution of each node on the path to the propagation of the event. The propagation contribution comprehensively considers factors such as the node's position in the path, edge weights, and the number of paths passing through that node. Different node positions in the path have different impacts on event propagation; edge weights reflect the tightness of the connection between nodes, with higher edge weights indicating stronger propagation effects. The more paths passing through a node, the more important that node is in the event propagation. Nodes are ranked according to their propagation contribution, and the node with the highest propagation contribution is most likely a candidate for the source location of the meteorological extreme event.

[0179] In one embodiment, step S550 may include the following steps S551 to S559:

[0180] Step S551: Perform topological sorting on the spatiotemporal propagation directed graph to generate a linear sequence of nodes, ensuring that all directed edges point from earlier nodes to later nodes in the sequence.

[0181] Topological sorting is an algorithm for sorting directed acyclic graphs (DAGs). It arranges the nodes in the graph into a linear sequence such that all directed edges point from earlier nodes to later nodes. Topological sorting of a spatiotemporally propagating directed graph generates a linear sequence of nodes, providing a sequential basis for subsequent dynamic programming calculations. Through topological sorting, it is ensured that when calculating the propagation contribution of nodes, the contribution of earlier nodes is calculated first, followed by the contribution of later nodes, conforming to the logical order of event propagation.

[0182] Step S552: Based on the topology sorting results and the weighted adjacency matrix, calculate the maximum weighted path from each node to all reachable nodes using dynamic programming, and record the node sequence, cumulative weight value, and path length in the path.

[0183] Dynamic programming is an algorithm that decomposes a complex problem into subproblems and uses the solutions to the subproblems to solve the original problem. Based on topological sorting and a weighted adjacency matrix, dynamic programming is used to calculate the maximum weighted path from each node to all reachable nodes. For each node, considering all its reachable subsequent nodes, the path weight from that node to its subsequent nodes is calculated based on the edge weights in the weighted adjacency matrix. By comparing the weights of different paths, the maximum weighted path is identified. The node sequence, cumulative weight value, and path length in the path are recorded; this information helps in subsequent analysis of the node's propagation contribution.

[0184] Step S553: ​​Construct a path-node association matrix. The rows of the matrix represent the propagation path, the columns represent the nodes, and the element values ​​are the position indices of the nodes in the corresponding paths. When a node is not in a path, the element value is -1.

[0185] A path-node association matrix is ​​used to represent the relationship between propagation paths and nodes. Rows in the matrix correspond to different propagation paths, and columns correspond to different nodes. Element values ​​are the position indices of nodes within their respective paths. If a node is in the path, the element value is its sequential number within the path (starting from 1); if a node is not in the path, the element value is -1. By constructing a path-node association matrix, it is easy to query the position information of each node in different propagation paths, providing fundamental data for calculating the propagation contribution of each node.

[0186] For each propagation path, traverse the nodes in the path and assign them position indices according to their order in the path. Create a path-node association matrix, where the number of rows equals the number of propagation paths and the number of columns equals the number of nodes. For each element in the matrix, assign a value based on whether the node is in the corresponding path and its position index within the path.

[0187] Step S554: Calculate the path coverage of each node, which is the ratio of the number of propagation paths containing that node to the total number of propagation paths, and generate a path coverage vector.

[0188] Path coverage reflects the degree to which a node is covered across all propagation paths. Path coverage for each node is calculated by dividing the number of propagation paths containing that node by the total number of propagation paths. The path coverage of each node is then combined to generate a path coverage vector. Higher path coverage generally indicates greater importance of the node in event propagation, as more propagation paths pass through it.

[0189] For each node, iterate through the path-node association matrix and count the number of elements in the column containing that node that are not -1. This count represents the number of propagation paths containing that node. Divide this number by the total number of propagation paths to obtain the path coverage for that node. Arrange the path coverage of each node in order of node sequence to generate a path coverage vector.

[0190] Step S555: Calculate the sum of the weights of the incoming edges and the sum of the weights of the outgoing edges for each node, and generate the in-degree centrality and out-degree centrality parameters.

[0191] The sum of incoming edge weights refers to the sum of the weights of all edges pointing to a node, while the sum of outgoing edge weights refers to the sum of the weights of all edges originating from a node. In-degree centrality and out-degree centrality parameters are defined based on the sum of incoming and outgoing edge weights, respectively. In-degree centrality reflects a node's ability to receive influence from other nodes, while out-degree centrality reflects a node's ability to propagate influence to other nodes. By calculating the sum of incoming and outgoing edge weights for each node, in-degree centrality and outgoing degree centrality parameters are generated. These two parameters help to comprehensively evaluate the role of a node in event propagation.

[0192] Step S556: Construct a node propagation influence diffusion model, taking path coverage, in-degree centrality, and out-degree centrality as inputs, and calculate the initial propagation contribution of the nodes using a preset weighted formula.

[0193] The node propagation influence diffusion model comprehensively considers factors such as path coverage, in-degree centrality, and out-degree centrality to calculate the initial propagation contribution of a node. A pre-defined weighting formula assigns different weights based on the importance of these factors to the node's propagation contribution. Path coverage, in-degree centrality, and out-degree centrality are used as inputs, and the initial propagation contribution of a node is calculated through a weighted summation. This model can more comprehensively evaluate the role of nodes in the propagation of extreme meteorological events.

[0194] Step S557: Perform path attenuation correction on the initial propagation contribution. Set the attenuation coefficient according to the position of the node in the propagation path. The attenuation coefficient of the starting node of the path is less than that of the ending node.

[0195] Path decay correction considers the impact of a node's position in the propagation path on its propagation contribution. A decay coefficient is set based on the node's position in the propagation path; the decay coefficient for nodes at the beginning of the path is smaller than that for nodes at the end, because nodes at the beginning of the path may play a more crucial role in the initiation and propagation of the event, and their propagation contribution should be relatively larger. Path decay correction is applied to the initial propagation contribution by multiplying it by the decay coefficient to obtain the corrected propagation contribution, making the calculation of the propagation contribution more reasonable.

[0196] Step S558: Multiply the initial propagation contribution by the attenuation coefficient to generate the final propagation contribution comprehensive score.

[0197] Multiplying the initial propagation contribution by the attenuation coefficient yields the final comprehensive propagation contribution score. This score comprehensively considers factors such as the node's position in the propagation path, path coverage, in-degree centrality, and out-degree centrality, thus more accurately reflecting the node's actual contribution to the propagation of extreme meteorological events. The final comprehensive propagation contribution score allows for a more reasonable ranking and comparison of nodes.

[0198] Step S559: Sort all nodes in descending order according to the comprehensive score of propagation contribution, and select the node with the highest score as the candidate source location of the propagation of the extreme meteorological event.

[0199] All nodes are sorted in descending order based on their final propagation contribution score. The node with the highest score represents the largest contribution to the propagation of the extreme weather event and is the most likely candidate for the source location of the event. This method allows for the selection of the most probable source node from a large pool of nodes. All nodes are then sorted from highest to lowest based on their final propagation contribution score. The node ranked first in this sorting is selected as the candidate for the source location of the extreme weather event.

[0200] Step S560: Perform multipath verification on the candidate propagation source location, calculate the probability confidence of the candidate location as the propagation source, and determine the candidate location as the final propagation source location of the extreme meteorological event when the confidence exceeds the preset confidence threshold.

[0201] Multipath validation is used to further confirm the reliability of candidate propagation source locations. The probability confidence level of a candidate location as a propagation source is calculated, considering factors such as its occurrence in multiple propagation paths and path weights. The probability confidence level reflects the likelihood of the candidate location being a propagation source. When the confidence level exceeds a preset threshold, it indicates that the candidate location has a high degree of credibility as a propagation source, and it is thus determined as the final propagation source location of the extreme meteorological event.

[0202] Step S570: Calculate the propagation influence parameters of the propagation source location. The propagation influence parameters include the propagation range radius, propagation speed coefficient, and influence attenuation index, and generate spatial propagation source identification results containing the propagation source coordinates and influence parameters.

[0203] The propagation radius refers to the maximum spatial range that the propagation source can influence, typically calculated as the distance to the farthest affected location centered on the propagation source. The propagation velocity coefficient reflects the speed at which extreme meteorological events propagate from the source to other locations; it can be calculated by analyzing the time difference and distance between nodes along the propagation path. The influence decay index describes the degree to which the influence of extreme meteorological events decreases with increasing distance from the propagation source. By calculating the propagation influence parameters at the source location and combining the source coordinates with these influence parameters, a spatial propagation source identification result is generated. This result comprehensively describes the propagation source and influence characteristics of extreme meteorological events.

[0204] Figure 2 A hardware entity diagram of a computer system provided as an embodiment of the present invention, such as... Figure 2 As shown, the hardware entity of the computer system 1000 includes a processor 1001 and a memory 1002, wherein the memory 1002 stores a computer program that can run on the processor 1001, and the processor 1001 executes the program to implement the steps in the method of any of the above embodiments.

Claims

1. A method for identifying spatial propagation sources based on extreme meteorological events, characterized in that, The method includes: Historical meteorological observation data of the study area are reconstructed into a time series to generate a meteorological data sequence set containing multiple meteorological indicators. The meteorological data sequence set has a uniform time sampling interval and spatial resolution. Based on a preset extreme event determination strategy, the meteorological data sequence set is processed to extract events, resulting in a set of meteorological extreme event sequences with spatiotemporal markers. This set includes the event occurrence time and corresponding spatial coordinates. Specifically, the process involves: inputting the meteorological data sequence set into an extreme event detection module; performing threshold crossing detection on each meteorological indicator sequence according to a preset multi-factor joint threshold rule to generate a preliminary extreme event candidate set, which includes meteorological data points exceeding the indicator threshold and their corresponding timestamps; performing spatial continuity analysis on the preliminary extreme event candidate set to identify clusters of extreme data points occurring simultaneously in adjacent spatial locations within a preset time window, generating event units with spatial clustering characteristics, where each event unit includes the spatial range covered by the cluster. The event unit is analyzed in several ways: 1) Calculate the spatial influence radius parameter based on the spatial range of the event unit, and 2) Calculate the event intensity index based on the duration interval. The event intensity index is positively correlated with the spatial influence radius parameter and the duration interval length. 3) Perform temporal continuity verification on the event unit, merge adjacent event units that are temporally continuous and spatially overlapping, and generate a set of meteorological extreme event sequences with unique event identifiers. 4) Add spatiotemporal marker information to each event unit in the set of meteorological extreme event sequences. The spatiotemporal marker information includes the event start time, end time, spatial center point coordinates, and event intensity index. 5) Establish an event attribute index table, and associate the spatiotemporal marker information with the unique identifier of the event unit for storage. The event attribute index table supports fast event retrieval operations based on spatiotemporal range. Spatiotemporal correlation strength is calculated on the set of meteorological extreme event sequences to generate a spatiotemporal correlation strength matrix describing the degree of correlation between events at different spatial locations. The element values ​​of the spatiotemporal correlation strength matrix represent the event synchronization level of the corresponding spatial location pair. The spatiotemporal correlation strength matrix is ​​subjected to a significance test, and the correlation strength elements that pass the significance threshold test are selected to construct a significant correlation set of meteorological extreme events; Based on the significant association set, reverse tracing analysis of the propagation path is performed to determine the spatial source location of the meteorological extreme event and the corresponding propagation influence parameters, and to generate spatial propagation source identification results containing the source coordinates and influence parameters. The step of performing threshold crossing detection on each meteorological indicator sequence according to a preset multi-element joint threshold rule to generate a preliminary extreme event candidate set includes: The time series curves of each meteorological indicator are extracted from the meteorological data sequence set. Each time series curve is decomposed by the ensemble empirical mode decomposition method to obtain a set of intrinsic mode functions containing different frequency components and trend components. Sliding window extremum detection is performed on each intrinsic modulus function and trend component to identify the local maxima and local minima of each component within the sliding window, generating a multi-component extremum set. The extreme points in the multi-component extreme point set are compared with the threshold range of the corresponding component, and abnormal extreme points that exceed the threshold range are filtered out. The timestamp, component type and corresponding meteorological index value of the abnormal extreme points are recorded. Consistency verification is performed on abnormal extreme points of different components at the same timestamp. When more than a preset number of components have abnormal extreme points at the same timestamp, the timestamp is marked as a joint abnormal moment. A time continuity analysis is performed on the joint anomaly moments, and consecutively occurring joint anomaly moments are merged into anomaly periods. The cumulative deviation of each meteorological indicator within the anomaly period is calculated. An abnormal event score is generated based on the duration and cumulative deviation of the abnormal period. Abnormal periods with scores exceeding the event score threshold are retained, and a preliminary extreme event candidate set containing the time interval and corresponding spatial location of the abnormal period is generated.

2. The method according to claim 1, characterized in that, The step of performing spatial continuity analysis on the preliminary extreme event candidate set to identify clusters of extreme data points that occur simultaneously in adjacent spatial locations within a preset time window includes: Construct a spatial grid index structure for the study area, map the abnormal data points in the preliminary extreme event candidate set to the corresponding grid cells, and generate a spatial grid event distribution with time stamps. Calculate the global spatial autocorrelation index of the spatial grid event distribution to determine whether there are significant clustering characteristics in the spatial distribution of extreme data points; When the global spatial autocorrelation index exceeds a preset clustering threshold, the spatial grid event distribution is clustered using a density-based spatial clustering algorithm. A distance threshold and a density threshold are set, and grid cells with a spatial distance less than the distance threshold and a number of data points exceeding the density threshold are aggregated into an initial extreme data point cluster. For each initial extreme data point cluster, boundary identification processing is performed. The spatial contour boundary of the cluster is calculated using the convex hull algorithm. The coordinates of the grid cells on the boundary are extracted as vertices of the cluster's spatial range to generate a boundary-optimized extreme data point cluster. Calculate the spatial centroid coordinates and coverage area parameters of the boundary optimization extreme data point cluster, and combine them with the temporal distribution characteristics of abnormal data points within the cluster to generate a characteristic extreme data point cluster containing start timestamp, end timestamp, spatial centroid coordinates, and coverage area parameters. Spatial overlap analysis is performed on the characteristic extreme data point clusters. When the spatial overlap area of ​​two characteristic extreme data point clusters exceeds the overlap threshold, they are determined to be spatially correlated clusters. Temporal correlation is then performed on the spatially correlated clusters by comparing the overlap interval length between the start and end timestamps of the clusters. When the overlap interval length exceeds the temporal correlation threshold, the two spatially correlated clusters are merged to generate a spatiotemporally correlated extreme data point cluster.

3. The method according to claim 1, characterized in that, The calculation of the spatiotemporal correlation strength of the meteorological extreme event sequence set to generate a spatiotemporal correlation strength matrix describing the degree of correlation between events at different spatial locations includes: The spatial location coordinates and time interval information of all event units are extracted from the set of meteorological extreme event sequences to construct an event spatiotemporal distribution dataset. The event spatiotemporal distribution dataset contains the event occurrence time series and the corresponding event intensity index for each spatial location. The event spatiotemporal distribution dataset is subjected to multi-scale spatial gridding processing, dividing the study area into multiple grid cell levels with different resolutions to form a multi-scale spatial analysis framework. Each grid cell level contains grid cells at the corresponding resolution and event statistical features within the cells. Under each scale grid of the multi-scale spatial analysis framework, the event time difference matrix between any two grid cells is calculated. The element values ​​of the event time difference matrix are statistical characteristic parameters of the absolute difference in the occurrence time of events in the corresponding grid cell pair. Based on the event time difference matrix and event intensity index at each scale, the spatiotemporal correlation coefficient is calculated at the corresponding scale. The spatiotemporal correlation coefficient is composed of a weighted combination of the synchronicity of event occurrence time and the correlation of event intensity. The fusion weights are set according to the event density of grid cells at each scale, and the spatiotemporal correlation coefficients at different scales are weighted and fused to generate cross-scale fused spatiotemporal correlation coefficients. The spatiotemporal correlation coefficients obtained from the cross-scale fusion are standardized to generate the final spatiotemporal correlation strength matrix. Each element in the matrix corresponds to the degree of correlation of a spatial location pair identified by both the row index and the column index.

4. The method according to claim 3, characterized in that, The process of performing multi-scale spatial gridding on the event spatiotemporal distribution dataset divides the study area into multiple grid cell levels of different resolutions, forming a multi-scale spatial analysis framework, including: Digital elevation model data of the study area are extracted from the geographic information database. A terrain complexity index is generated by calculating the elevation standard deviation and the mean slope. The terrain complexity index is positively correlated with the elevation standard deviation and the mean slope. The basic grid resolution parameters are determined based on the terrain complexity index. The higher the terrain complexity index, the higher the basic grid resolution parameters. Based on the basic grid resolution parameters, multiple resolution levels are generated according to the hierarchical progression to form a resolution sequence. The resolution of adjacent levels increases by a preset ratio, so that the grid units of each level maintain a nested relationship in terms of spatial coverage. The study area is meshed layer by layer according to the resolution sequence to generate mesh unit levels of different resolutions. Each mesh unit level contains a set of mesh units and unit boundary coordinates at the corresponding resolution. Each grid cell level is spatially encoded, and a unique spatial index code is assigned to each grid cell through a hierarchical spatial index. The spatial index code includes the level identifier and the two-dimensional coordinate offset of the cell within the level. The event units in the event spatiotemporal distribution dataset are mapped to each grid unit level according to their spatial location. The distribution characteristics of event occurrence frequency, average event intensity index and event duration in each grid unit are statistically analyzed to generate a grid event statistical feature table containing spatial index codes.

5. The method according to claim 4, characterized in that, The calculation of the event time difference matrix between any two grid cells within each scale of the multi-scale spatial analysis framework includes: Extract the event occurrence time series of the target scale grid cell level in the multi-scale spatial analysis framework from the grid event statistical feature table. The event occurrence time series includes the event start timestamp and corresponding event intensity index of all grid cells at that level. The event occurrence time series of each grid cell is aligned with the time series, and the event state values ​​of missing time points are supplemented by time interpolation method to generate an event state series with a uniform time interval. The time interval of the event state series is consistent with the event sampling frequency of the grid at that scale. Calculate the cross-correlation function of event state sequences of any two grid cells, identify the time lag corresponding to the maximum value of the cross-correlation function, and use the time lag as the characteristic time difference parameter of the two grid cells. The characteristic time difference parameter characterizes the propagation time feature of the event between the two grid cells. An initial time difference matrix is ​​constructed based on the characteristic time difference parameters of all grid cell pairs. The rows and columns of the matrix correspond to different grid cells, and the element values ​​are the characteristic time difference parameters of the corresponding grid cell pairs. The initial time difference matrix is ​​normalized by time scale, and the matrix element values ​​are converted into proportions relative to the average duration of events under that scale grid to generate a normalized time difference matrix. The spatial autocorrelation index of the normalized time difference matrix is ​​calculated. When the spatial autocorrelation index is lower than the preset correlation threshold, the matrix is ​​spatially smoothed. Local noise interference is eliminated by the neighborhood weighted average method to generate event time difference matrices at each scale.

6. The method according to claim 1, characterized in that, The step involves performing a significance test on the spatiotemporal correlation strength matrix, selecting correlation strength elements that pass the significance threshold test, and constructing a significant correlation set for meteorological extreme events, including: A random alternative dataset with the same temporal and spatial distribution characteristics as the set of meteorological extreme event sequences is constructed. The random alternative dataset is generated by randomly rearranging the original event time series. Based on the random substitution dataset, the spatiotemporal correlation strength calculation process is repeated to generate multiple random correlation strength matrices. The number of random correlation strength matrices is determined by a preset number of repeated samplings. From each random association strength matrix, the association strength value corresponding to the position of each element in the original spatiotemporal association strength matrix is ​​extracted, and a random distribution model of each element is constructed. The random distribution model includes the probability density function and the cumulative distribution function of the association strength value. Calculate the quantile of each element value in the original spatiotemporal correlation strength matrix in the corresponding random distribution model, and generate a significance p-value, which represents the probability that the element value is generated by a random process; Filter out the association strength elements with a significance p-value less than a preset significance threshold, and record the corresponding grid cell pairs and association strength values ​​of these elements; Spatial connectivity analysis was performed on the selected correlation strength elements to construct a correlation network. The set of correlation strength elements that form connected components was retained, and isolated correlation strength elements were deleted to construct the final significant correlation set of meteorological extreme events.

7. The method according to claim 6, characterized in that, The construction of a random alternative dataset with the same temporal and spatial distribution characteristics as the set of meteorological extreme event sequences includes: The time interval information and spatial location coordinates of all event units are extracted from the set of meteorological extreme event sequences. The time periodic characteristics of the event occurrence are calculated by the periodogram analysis method, and the spatial density distribution characteristics of the event occurrence are calculated by the kernel density estimation method. An event time point process model is constructed based on the aforementioned time period characteristics, and the intensity function of the event time point process model is consistent with the time distribution characteristics of the original event. An event space sampling model is constructed based on the aforementioned spatial density distribution characteristics to generate random spatial coordinates that maintain the original spatial density distribution. The event sequence at each spatial location is time-rearranged, keeping the number of events unchanged but redistributing the event occurrence times according to the event time point process model, to generate a time-rearranged event sequence; Spatial location perturbation processing is applied to the time-rearranged event sequence by adding random perturbation terms that follow a preset distribution to the original spatial location coordinates to generate a spatial perturbation event sequence. The spatial perturbation event sequence is subjected to spatiotemporal consistency verification to check whether the autocorrelation characteristics of the time series are consistent with the original sequence. The spatial distribution characteristics are checked by the spatial autocorrelation index. Event units that do not meet the consistency requirements are deleted, and the final random alternative dataset is generated.

8. A computer system comprising a memory and a processor, the memory storing a computer program executable on the processor, characterized in that, When the processor executes the program, it implements the steps of the method according to any one of claims 1 to 7.

Images