A Regional Cooperative Storm Surge Prediction Method Using Coupled Graph Attention and Gated Recurrent Networks

Abstract

This invention relates to the field of storm surge disaster monitoring and intelligent forecasting technology, and provides a regional collaborative storm surge rise forecasting method using coupled graph attention and gated cyclic networks. The method includes: using a hybrid wind field-driven two-dimensional hydrodynamic model formed by fusing typhoon parameterization and reanalysis data to obtain storm surge rise residual field data; registering the rise residual field data with the observation time series in time and space, and using it as a physical prior input; constructing a dynamic edge-weighted graph structure based on a spatiotemporal causal correlation analysis method to adaptively update the dependencies between multiple observation stations; using the physical prior and key features as node inputs, fusing the graph attention mechanism of GATv2 and the temporal feature extraction capability of GRU to perform spatiotemporal coupled modeling of storm surge rise; conducting multi-site joint prediction through regional collaboration, performing model evaluation and result optimization, and outputting storm surge rise forecast sequences with a 15-minute time step and a lead time of 15 minutes to 4 hours.

Description

Technical Field

[0001] This invention belongs to the field of storm surge disaster monitoring and intelligent forecasting technology, specifically involving a regional collaborative storm surge rise forecasting method using coupled graph attention and gated cyclic networks. Background Technology

[0002] Storm surges are sudden coastal disasters triggered by tropical or extratropical cyclones. Their core hazard lies in the superposition of astronomical tides and storm surges, resulting in water levels exceeding warning levels, which can cause seawall breaches, urban flooding, and ecosystem damage. Existing storm surge prediction methods mainly fall into three categories: numerical simulation, traditional deep learning, and graph neural networks. While numerical models can characterize the physical processes, they are computationally intensive and time-consuming, making them unsuitable for rapid early warning. Traditional deep learning methods, although computationally efficient, can only capture the temporal characteristics of a single station, ignoring the spatial propagation relationships of storm surges across multiple stations. Existing graph neural networks attempt to characterize spatial dependencies through graph structures, but they often use static edge weights, lacking physical constraints and dynamic adaptability, leading to increased prediction errors during typhoon peak periods. In summary, existing methods struggle to achieve high-precision, long-term, and robust predictions of storm surge increases at multiple stations in semi-enclosed sea areas.

[0003] To address the aforementioned issues, this invention proposes a regional collaborative storm surge prediction method that couples graph attention with a gated cyclic network. This method integrates physics-driven and data-driven mechanisms, using MIKE21 HDFM numerical simulation results as hydrodynamic prior constraints to construct a graph attention-gated cyclic model capable of adaptively updating edge weights. This enables dynamic spatial dependency modeling and deep fusion of spatiotemporal features among multiple storm surge observation points. Summary of the Invention

[0004] This invention aims to overcome the problems of long computation time, insufficient spatial dependence characterization, and limited prediction accuracy in existing storm surge prediction methods. It proposes a regional collaborative storm surge prediction method using coupled graph attention and gated recurrent networks. By combining physical numerical simulation with data-driven learning and introducing a dynamic graph attention mechanism and temporal memory structure, a spatiotemporal collaborative prediction framework that reflects the dynamic characteristics of typhoons is constructed.

[0005] The technical solution adopted by this invention to solve the above-mentioned technical problems is: a regional collaborative storm surge prediction method using coupled graph attention and gated recurrent networks, comprising:

[0006] Based on the hybrid wind field driven two-dimensional hydrodynamic model formed by the fusion of typhoon parameterization and reanalysis data, storm surge water increase residual field data is obtained. The water increase residual field data is registered with the observation time series in time and space and used as physical prior input.

[0007] A dynamic edge-weighted graph structure is constructed based on the spatiotemporal causal correlation analysis method to adaptively update the dependencies between multiple observation stations;

[0008] By using physical priors and key features as node inputs, and integrating the graph attention mechanism of GATv2 with the temporal feature extraction capability of GRU, spatiotemporal coupled modeling of storm surge rise is performed.

[0009] By conducting joint forecasts at multiple sites through regional collaboration, model evaluation and result optimization are performed, and storm surge forecast sequences with a 15-minute time step and a lead time of 15 minutes to 4 hours are output.

[0010] As a preferred embodiment, the specific steps for obtaining storm surge residual field data are as follows:

[0011] A hybrid wind field is constructed by integrating the typhoon parametric wind field model with reanalysis data. First, the typhoon parametric wind field is calculated. Based on the wind field parametric model, the wind pressure structure of the typhoon core is determined, as shown in the following formula:

[0012] ,

[0013] in, Distance from the typhoon center The air pressure at that location The radius of the distance from the center of the typhoon. The central pressure of the typhoon. For ambient air pressure, The radius of maximum wind speed. Here, exp is the shape parameter, and exp is the exponential function;

[0014] Calculation of typhoon wind speed based on gradient wind formula:

[0015] ,

[0016] in, Distance from the typhoon center Gradient wind speed at the location, air density, For Coriolis parameters;

[0017] Furthermore, global reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) were selected, with a spatial resolution of 0.25°×0.25° and a temporal resolution of 1 hour. The latitude and longitude range of the study area was cropped, and 10m wind speed, wind direction, and sea level pressure data were extracted as a large-scale background circulation field.

[0018] A weighted smoothing algorithm is used to fuse the two data in the typhoon's core region. The weight of the typhoon parameterized wind field is set at 0.8-1.0, and the weight of the wind field in the reanalysis data is set at 0.0-0.2; in the typhoon field... The typhoon parameterized wind field weights are set to 0.0-0.5, and the wind field weights of the reanalysis data are set to 0.5-1.0. After fusion, a 5×5 Gaussian filter is used to eliminate spatial heterogeneity and generate hourly mixed wind fields with a resolution of 0.1°.

[0019] As a preferred embodiment, the specific steps for registering the water-enhanced residual field with the observation time series in time and space are as follows:

[0020] After converting all data to UTC, the straight-line distance between the tide gauge station and the typhoon center is calculated using the Haversine function, based on the latitude and longitude of the tide gauge station and the typhoon center. The formula is as follows:

[0021] ,

[0022] Where R is the Earth's radius. This represents the latitude difference between the station and the typhoon center. This is the difference in longitude between the station and the typhoon center. , These are the latitudes of the tide gauge station and the typhoon center, respectively.

[0023] Furthermore, piecewise constants are used to smooth typhoon level data, and missing time points are filled with prior raw values. For the maximum wind speed near the center and the central pressure, dynamic window Savitzky-Golay filtering combined with piecewise linear interpolation is used for time-series preprocessing, and the filter window length is set. ,in The length is half the window length, dynamically set to 1-3. For odd numbers between 3 and 7, the fitting polynomial has an order of 2, and the filtering formula is as follows:

[0024] ,

[0025] Where t is the current time, and i is the number of deviations at the current time t. These are the filter coefficients. The original data within the window. This is the filtered data;

[0026] For the filtered 1-hour data, piecewise linear interpolation is used to generate 15-minute data. Time data is , Time data is ,but The interpolation result at time t is: Cubic spline interpolation is used to improve the resolution of the water level rise residual sequence within the interval. Let m be the number of adjacent intervals in the sequence, and construct a cubic polynomial. Apply interpolation conditions , Smooth conditions , and boundary conditions , The coefficients are obtained by solving a system of linear equations. , , , This leads to the generation of a 15-minute resolution water replenishment residual sequence.

[0027] As a preferred implementation, the specific steps for constructing a dynamic edge-weighted graph structure based on the spatiotemporal causal correlation analysis method are as follows:

[0028] Based on the Liang-Kleeman information flow, a dynamic weighted graph structure is constructed. First, a graph G=(V,E,W) is defined, where V={A,B,C,D,E} are the tide gauge station nodes, E are the causal edges between nodes, and W is the edge weight matrix.

[0029] Secondly, the information flow intensity between stations is calculated, and the causal impact is quantified using the Liang-Kleeman information flow theory for the time series of water increase residuals at each station, as shown in the following formula:

[0030] ,

[0031] Where A is the driving variable in the time series, and B is the response variable in the time series. Let A be the time derivative. For covariance, For variance;

[0032] Introducing a time lag parameter τ, we calculate the information flow intensity under different τ values ​​and select one that makes the information flow intensity more stable. The largest τ is used as the optimal lag, and causal relationships are filtered through information flow values;

[0033] A graph neural network based on GATv2 is constructed. Features are input to node s and its neighboring nodes g. Attention coefficients are calculated using a linear transformation and the LeakyReLU activation function, as shown in the following formula:

[0034] ,

[0035] ,

[0036] in, Scoring unnormalized attention Let T be the attention vector, and T be the matrix transpose symbol. For trainable weight matrix, , Let be the feature vectors of nodes s and g. Let be the neighborhood of node s. is the attention coefficient, softmax is the normalization function, exp is the exponential function, and l represents the loop index.

[0037] As a preferred implementation, the specific steps of fusing the graph attention mechanism of GATv2 with the temporal feature extraction capability of GRU are as follows:

[0038] The features of neighboring nodes are weighted and aggregated using attention coefficients and activation functions are applied. The multi-head attention outputs are averaged among the heads to generate a node space representation.

[0039] A GRU sequence model is constructed, defining reset and update gates. Based on the current input and the hidden state from the previous time step, the gate value and candidate state are obtained, and element-wise gating state updates are performed. The hidden state is calculated, and the hidden state of the last GRU layer is used as the input linear layer to map the predicted value. The storm surge prediction value is output for each station, as shown in the following formula:

[0040] ,

[0041] in This is the output layer weight matrix. For output layer bias terms, Let T be the hidden state vector at time step T. This is a predicted value;

[0042] The model was trained by setting parameters, using mean squared error as the loss function, and updating the model parameters using the Adam optimizer. The model overfitting was evaluated using the validation set, and the number of GATv2 attention heads and GRU layers was adjusted accordingly. The trained model was then applied to an independent test set to generate preliminary storm surge prediction results.

[0043] As a preferred implementation method, the specific steps for model evaluation and result optimization are as follows:

[0044] The root mean square error, mean absolute error, and coefficient of determination are calculated to evaluate the consistency between the predicted and measured sequences. Under the same data and evaluation process, the differences of the above indicators between the TGAT model and SVR, LSTM, and atLSTM are compared to obtain the horizontal comparison results of different models. The distribution characteristics in the spatial and temporal dimensions are analyzed by statistically and visually predicting residuals, focusing on locating high error intervals and checking the quality of input data. Local regression smoothing technology is used to correct the temporal neighborhood of samples identified as having high errors.

[0045] Beneficial effects

[0046] Compared with the prior art, the present invention has the following advantages:

[0047] 1. Integrating the advantages of both physical and data-driven approaches, and using numerical simulation results as prior constraints, it improves the prediction stability under extreme weather conditions, enables dynamic spatial dependency modeling of multiple sites, can adaptively capture changes in node correlation intensity caused by wind field and air pressure variations, takes into account both short-term and long-term time series memory, improves the accuracy of multi-time-effect water increase forecasts, supports regional collaborative early warning and emergency decision-making, and has good generalization and engineering usability.

[0048] 2. Breaking through the limitations of traditional static graph structures, this model dynamically adjusts edge weights based on information transmission intensity and typhoon development stages to achieve real-time updates of spatial dependencies and physical rationality constraints. It utilizes GATv2's multi-head attention mechanism to distinguish the impact priorities of different observation stations and combines GRU's gated memory structure to suppress the accumulation of temporal errors, enabling collaborative learning of short-term and long-term features. The 4-hour prediction accuracy can be maintained above 0.91. Numerical simulation data ensures that the input conforms to hydrodynamic laws. Physically impossible paths are excluded during the graph structure construction stage. Typhoon feature parameters are introduced during model training to enhance physical interpretability, significantly improving the model's stability and generalization performance.

[0049] 3. By combining physical priors, dynamic graph structures, and temporal memory, the traditional numerical model overcomes the problems of long computation time and insufficient expressive power of static graph networks. It significantly improves the accuracy, stability, and real-time performance of storm surge forecasting in semi-enclosed shallow sea areas. It can be directly applied to coastal disaster prevention and mitigation, seawall reinforcement planning, and emergency dispatch decision-making, providing efficient and reliable technical support for intelligent early warning of storm surge disasters. Attached Figure Description

[0050] Figure 1 This is a flowchart of the present invention;

[0051] Figure 2 This is the location of the research area in this invention;

[0052] Figure 3 This is a comparison chart of the error distribution of the present invention, wherein (a) is the overall error histogram and (b) is the error violin plot of each time step and station;

[0053] Figure 4 This invention provides a multi-station comparison of the prediction performance of TGAT and GRU within different prediction time ranges.

[0054] Figure 5 This is a performance comparison of different graph structures of the present invention, wherein 1 is a distance static graph, 2 is a similarity static graph, 3 is an information flow static graph, and 4 is a TGAT dynamic graph;

[0055] Figure 6 This invention provides a heatmap of the attention weights of sites at different stages of a typhoon. Detailed Implementation

[0056] To make the technical means, creative features, achieved objectives, and effects of this invention easier to understand, the invention is further described below with reference to specific embodiments. However, the following embodiments are merely preferred embodiments of this invention and not all of them. Other embodiments obtained by those skilled in the art based on the embodiments described herein without creative effort are all within the protection scope of this invention.

[0057] To achieve the above objectives, such as Figure 1 As shown, a regional collaborative storm surge forecasting method using coupled graph attention and gated recurrent networks is presented, the method comprising the following steps:

[0058] Step S1: Data acquisition and numerical simulation of storm surge residuals;

[0059] Step S2: Data preprocessing and feature selection;

[0060] Step S3: TGAT model construction and training;

[0061] Step S4: Model evaluation and result optimization;

[0062] Step S5: Output and visualization of results.

[0063] This method is implemented in the order of S1–S5, and its overall process is as follows:

[0064] Based on the hybrid wind field driven two-dimensional hydrodynamic model formed by the fusion of typhoon parameterization and reanalysis data, storm surge water increase residual field data is obtained. The water increase residual field data is registered with the observation time series in time and space and used as physical prior input.

[0065] A dynamic edge-weighted graph structure is constructed based on the spatiotemporal causal correlation analysis method to adaptively update the dependencies between multiple observation stations;

[0066] By using physical priors and key features as node inputs, and integrating the graph attention mechanism of GATv2 with the temporal feature extraction capability of GRU, spatiotemporal coupled modeling of storm surge rise is performed.

[0067] By conducting joint forecasts at multiple sites through regional collaboration, model evaluation and result optimization are performed, and storm surge forecast sequences with a 15-minute time step and a lead time of 15 minutes to 4 hours are output.

[0068] Combination Figure 2 – Figure 6 The specific implementation methods and effects of steps S1–S5 are explained respectively, with the specific steps of data acquisition and numerical simulation of storm surge residuals as follows:

[0069] A hybrid wind field is constructed by integrating the Holland typhoon model and ERA5. First, the Holland typhoon wind field is calculated. Based on the Holland parametric model, the wind pressure structure of the typhoon core is determined, as shown in the following formula:

[0070] ,

[0071] in, Distance from the typhoon center The air pressure at that location The radius of the distance from the center of the typhoon. The central pressure of the typhoon. The ambient air pressure is taken as 1010 hPa. The radius of maximum wind speed. Here, is the shape parameter, which can be 2.0, and exp is an exponential function;

[0072] Calculation of typhoon wind speed based on gradient wind formula:

[0073] ,

[0074] in, Distance from the typhoon center Gradient wind speed at the location, For air density, take 1.225 kg / m³. 3 , For Coriolis parameters;

[0075] Furthermore, the ERA5 global reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) was selected, with a spatial resolution of 0.25°×0.25° and a temporal resolution of 1 hour. The latitude and longitude range of the study area was cropped to 109.5°E–111.0°E and 19.5°N–20.5°N. 10m wind speed, wind direction and sea level pressure data were extracted as large-scale background circulation field.

[0076] Finally, a weighted smoothing algorithm is used to fuse the two in the typhoon's core region. The weighting for the Holland wind field is 0.8-1.0, and the weighting for the ERA5 wind field is 0.0-0.2; in the typhoon field... The Holland wind field weight is set to 0.0-0.5, and the ERA5 wind field weight is set to 0.5-1.0. After fusion, a 5×5 Gaussian filter is used to eliminate spatial heterogeneity and generate hourly mixed wind fields with a resolution of 0.1°.

[0077] The MIKE21 HDFM model parameters were configured, using an unstructured triangular mesh covering the area from 109°E to 111.5°E and from 19°N to 20.8°N. In the nearshore area, the mesh resolution was increased to 50–200 m to ensure the capture of shoreline irregularities and shoal topography. In the offshore area, the mesh resolution was gradually increased from 200 m to 1 km to balance computational accuracy and efficiency. The water depth data used was the EMODnet 2023 version 500 m resolution topographic dataset, which was interpolated and adapted to the model mesh. The average water depth in the Qiongzhou Strait area was set to 44 m, and the maximum water depth was set to 114 m.

[0078] The open boundary adopts the astronomical tidal boundary conditions provided by the TPXO9 global tidal model, including eight major tidal constituents: M2, S2, N2, K2, K1, O1, P1, and Q1. The surface drive uses the fused mixed wind field and ERA5 pressure field, and the bottom friction uses the Manning coefficient, which is 0.025 s / m for nearshore shoals. 1 / 3 0.018s / m at offshore location 1 / 3 The side boundaries adopt a no-slip condition;

[0079] Seven historical typhoon events from 2005 to 2024 were selected. The simulation period covered the critical period of typhoon impact, namely 24 hours before the typhoon center entered the study area and 24 hours after it left. The model time step was set to 60 seconds, the output frequency was 1 hour, and the output variables included total water level, flow velocity, and water increase residual.

[0080] The hourly water level rise residual sequences of five tide gauge stations were extracted from the model output. The coordinates of the tide gauge stations were determined based on the "Compilation of Historical Data of Chinese Ocean Stations". Outliers were removed using the 3σ criterion, and the mean water level rise residual of each station was calculated. with standard deviation Remove those smaller than or greater than For the sample, the two-station comparison method was used to verify consistency. For geographically adjacent stations, the correlation coefficient of their water increase residuals must be ≥0.85; otherwise, the model parameters were re-examined and the simulation was reworked.

[0081] Generate a user folder with the naming format "Typhoon_[Year]_[Typhoon Name]" to store mixed wind field data, MIKE model parameter files, and water rise residual sequences. The CSV file contains fields such as time, site ID, water rise residual, wind speed, and central pressure.

[0082] The specific steps of data preprocessing and feature selection are as follows:

[0083] First, all data is converted to UTC. Then, based on the latitude and longitude of the tide gauge station and the typhoon center, the straight-line distance between the tide gauge station and the typhoon center is calculated using the Haversine function, as shown in the following formula:

[0084] ,

[0085] Where R is the Earth's radius. This represents the latitude difference between the station and the typhoon center. This is the difference in longitude between the station and the typhoon center. , These are the latitudes of the tide gauge station and the typhoon center, respectively.

[0086] Furthermore, piecewise constants are used to smooth typhoon level data, and the most recent preceding raw values ​​are used to fill in missing moments. For the maximum wind speed near the center and the central pressure, dynamic window Savitzky-Golay filtering combined with piecewise linear interpolation is used for time series preprocessing, and the window length of the filter is set. ,in The length is half the window length, dynamically set to 1-3. For odd numbers between 3 and 7, the fitting polynomial has an order of 2, and the filtering formula is as follows:

[0087] ,

[0088] Where t is the current time, and i is the number of deviations at the current time t. These are the filter coefficients. The original data within the window. This is the filtered data;

[0089] For the filtered 1-hour data, piecewise linear interpolation is used to generate 15-minute data. Time data is , Time data is ,but The interpolation result at time t is: Cubic spline interpolation is used to improve the resolution of the water level rise residual sequence within the interval. Let m be the number of adjacent intervals in the sequence, and construct a cubic polynomial. Apply interpolation conditions , Smooth conditions , and boundary conditions , The coefficients are obtained by solving a system of linear equations. , , , This leads to the generation of a 15-minute resolution water replenishment residual sequence;

[0090] The distance between the tide gauge station and the typhoon center was interpolated using a fixed-window Savitzky-Golay filter and piecewise linear interpolation. The window length was fixed at 5, the polynomial order was 2, and the interpolation was performed using the same water level residual sequence method.

[0091] Typhoon level, maximum wind speed near the center, typhoon center pressure, distance between tide gauge station and typhoon center, and storm surge increase in the preceding time series were selected as key input features. The min-max method was used to normalize the feature data to the [0,1] interval, and the feature data were divided into training set, validation set and test set in a ratio of 7:1:2.

[0092] The specific steps for building and training the TGAT model are as follows:

[0093] like Figure 2 Taking the storm surge forecast of five tide gauge stations in the Qiongzhou Strait region (A, B, C, D, and E) as an example, a dynamic weighted graph structure is constructed based on Liang-Kleeman information flow. First, the graph G=(V,E,W) is defined, where V={A,B,C,D,E} are the nodes of the five tide gauge stations, E is the causal edge between nodes, and W is the edge weight matrix.

[0094] Secondly, the information flow intensity between stations is calculated, and the causal impact is quantified using the Liang-Kleeman information flow theory for the time series of water increase residuals at each station, as shown in the following formula:

[0095] ,

[0096] Where A is the driving variable in the time series, and B is the response variable in the time series. Let A be the time derivative. For covariance, For variance;

[0097] Taking storm surge forecasting in the Qiongzhou Strait region as an example, under the same evaluation process, such as Figure 5 As shown, four types of graph structures are constructed: 1. Distance static graph, 2. Similarity static graph, 3. Information flow static graph, and 4. TGAT dynamic graph. A comprehensive evaluation and horizontal comparison of the prediction results of the TGAT dynamic graph on the test set are conducted. The results show that the peak bias is effectively suppressed after adopting dynamic edge weights. On the test set, the TGAT model's RMSE for the 24-hour storm surge increase of Typhoon Kalmegi is 15.3 cm, MAE is 11.2 cm, and R... 2With a value of 0.89, it is superior to traditional models, which makes it easier to capture the dynamic causal relationship between tide gauge stations under different typhoon events. It is beneficial to overcome the limitations of traditional static graph structures and improve prediction accuracy. By capturing spatial features through the GATv2 module and temporal features through the GRU module, it achieves effective fusion of spatial-temporal multi-scale features and solves the technical problem that traditional methods cannot handle spatial correlation and temporal dependence at the same time.

[0098] Introducing a time lag parameter τ, we calculate the information flow intensity under different τ values, such as a range of -3 hours to +3 hours, with a step size of 15 minutes to set τ, and then selecting the appropriate setting. The largest τ is used as the optimal lag, and causal relationships are filtered by maximizing the information flow value;

[0099] Based on physical laws, impossible paths are eliminated, and an adaptive threshold is used to retain stable edges, such as the 20th percentile of non-zero information flow. The average information flow intensity is calculated as the edge weight, and a weight matrix W is generated.

[0100] Construct a graph neural network based on GATv2, setting two graph attention hidden layers with 16 neurons in each layer, and using two parallel attention heads in each layer;

[0101] For node s and its neighboring node g as input features, the attention coefficients are calculated using a linear transformation and the LeakyReLU activation function, as shown in the following formula:

[0102] ,

[0103] ,

[0104] in, Scoring unnormalized attention Let T be the attention vector, and T be the matrix transpose symbol. For trainable weight matrix, , Let be the feature vectors of nodes s and g. Let be the neighborhood of node s. Here, is the attention coefficient, softmax is the normalization function, exp is the exponential function, and l represents the loop index;

[0105] The features of neighboring nodes are weighted and aggregated using attention coefficients and activation functions are applied. The multi-head attention outputs are averaged among the heads to generate a node space representation.

[0106] Construct a GRU sequence model with 3 hidden layers, each with 16 neurons, and use the activation function as follows: Define reset and update gates. Based on the current input and the hidden state from the previous time step, obtain the gate value and candidate state, and perform element-wise gating state updates. Calculate the hidden state, take the hidden state of the last GRU layer, and map it to the linear layer as the predicted value. Output the storm surge prediction value for each station, as shown in the following formula:

[0107] ,

[0108] in This is the output layer weight matrix. For output layer bias terms, Let T be the hidden state vector at time step T. This is a predicted value;

[0109] The model was trained with the following parameters: learning rate 0.01, batch size 32, number of iterations 100, early stopping patience value 20, and early stopping threshold 0.0002. Mean squared error was used as the loss function, and the Adam optimizer was used to update the model parameters. The overfitting of the model was evaluated through the validation set, and the number of GATv2 attention heads and GRU layers was adjusted. The trained model was then applied to an independent test set to generate preliminary storm surge prediction results.

[0110] In this embodiment, the method first comprehensively evaluates the prediction results of the TGAT model on the test set; for the 2023 Typhoon Talim event, the average RMSE calculated for each tide gauge station is 15.3 cm, MAE is 11.2 cm, and R... 2 The value of 0.89 indicates that the TGAT model has high prediction accuracy; for example Figure 4 As shown, compared with GRU, TGAT has lower MAE and maximum error at all lead times, and the error hot zone converges, demonstrating the suppression effect of dynamic graph structure and regional collaborative modeling on peak stage error and long-term stability.

[0111] The specific steps for model evaluation and result optimization are as follows:

[0112] The root mean square error, mean absolute error, and coefficient of determination are calculated to evaluate the consistency between the predicted and measured sequences. Under the same data and evaluation process, the differences of the above indicators between the TGAT model and SVR, LSTM, and atLSTM are compared to obtain the horizontal comparison results of different models. The distribution characteristics in the spatial and temporal dimensions are analyzed by statistically and visually predicting residuals. The high error intervals, such as the peak period of typhoon passage, are located in the key area. The quality of the input data is checked in conjunction with the analysis. The temporal neighborhood of the samples identified as high error samples is corrected by local regression smoothing technology.

[0113] Local regression smoothing techniques are used to correct the predicted values ​​in these high-error regions, such as... Figure 3As shown, the corrected RMSE further decreased to 13.8 cm, MAE decreased to 9.9 cm, and R 2 The value was increased to 0.91, which further optimized the performance of the TGAT model, making it easier to comprehensively and objectively evaluate the performance of the TGAT model and overcoming the limitations of a single indicator.

[0114] The specific steps for outputting and visualizing the results are as follows:

[0115] By exporting the prediction results file and simultaneously outputting the model parameter file, including GATv2 weights and GRU weights, as well as the accuracy evaluation report, and by generating a comparison chart of water increase time series at multiple sites and a heat map of attention weights at different stages of typhoon passage, a visual result that intuitively presents the model’s temporal fit and spatial relationship evolution is obtained.

[0116] like Figure 6 As shown, the attention weight heatmaps for different stages of typhoon passage—6 hours before landfall, at landfall time, and 6 hours after landfall—use color intensity to represent the correlation strength between tide gauge stations, revealing the dynamic process and main path of tidal wave propagation. This helps to formulate differentiated early warning strategies and improve the scientific nature and pertinence of decision-making. The attention weight heatmaps for different stages of typhoon passage innovatively demonstrate the dynamic process of storm surge propagation, making complex physical mechanisms intuitive and easy to understand.

[0117] The embodiments of the present invention described above are subject to modification and change of method by those skilled in the art without departing from the embodiments and broader aspects of the present invention. The appended claims are intended to include all such modifications and changes of method that do not depart from the present invention.

Claims

1. A regional storm surge forecasting method coupling graph attention with gated recurrent networks, characterized in that, The method comprises the following steps: a two-dimensional hydrodynamic model driven by a mixed wind field formed by fusing typhoon parameterization and reanalysis data is used to obtain a storm surge residual field data, the storm surge residual field data is matched in time and space with observation time series, and the storm surge residual field data and the observation time series are used as physical priori input; a dynamic edge weight graph structure is constructed based on a spatiotemporal causal correlation analysis method, and a self-adaptive update is performed on the dependence relationship between multiple observation stations; a GATv2 graph attention mechanism and a GRU temporal feature extraction capability are fused, and a spatiotemporal coupling model is established for storm surge by taking the physical priori and key features as node input; a multi-station joint prediction is performed through regional cooperation, model evaluation and result optimization are performed, and a storm surge prediction sequence with a 15-minute time step and a 15-minute to 4-hour lead time is output.

2. The regional storm surge forecasting method of claim 1, wherein: The specific steps for obtaining the storm surge residual field data are as follows: a mixed wind field is constructed by fusing a typhoon parameterization wind field model and a reanalysis data wind field, the typhoon parameterization wind field is first calculated, the wind pressure structure of the typhoon core is determined based on the wind field parameterization model, and the formula is as follows: ， wherein, P is the atmospheric pressure at the center of the typhoon, R is the radius from the center of the typhoon, P0 is the atmospheric pressure at the center of the typhoon, P0 is the atmospheric pressure at the center of the typhoon, P0 is the atmospheric pressure at the center of the typhoon, Rmax is the radius of the maximum wind speed, A is a shape parameter, and exp is an exponential function. the typhoon wind speed is calculated based on the gradient wind formula: ， wherein, is the gradient wind speed at a distance of from the center of the typhoon, is the air density, is the Coriolis parameter; further, global reanalysis data with a spatial resolution of 0.25°x0.25° and a temporal resolution of 1 hour is selected, the latitude and longitude range of the research area is cropped, and 10m wind speed, wind direction and sea level pressure data are extracted as large-scale background circulation field; The two are fused by using weight smoothing algorithm. In the core area of typhoon , the weight of typhoon parameterization wind field is 0.8-1.0, and the weight of reanalysis data wind field is 0.0-0.2; in the outer field of typhoon , the weight of typhoon parameterization wind field is 0.0-0.5, and the weight of reanalysis data wind field is 0.5-1.0, after fusion, 5x5 Gaussian filter is used to eliminate spatial heterogeneity, and the mixed wind field with 0.1° resolution is generated every hour.

3. The regional storm surge forecasting method of claim 1, wherein: The specific steps for matching the storm surge residual field and the observation time series in time and space are as follows: all data are converted to coordinated universal time, the haversine function is used to calculate the straight-line distance between the tide gauge station and the typhoon center based on the latitude and longitude of the tide gauge station and the latitude and longitude of the typhoon center, and the formula is as follows: ， where R is the earth radius, is the latitude difference between the station and the typhoon center, is the longitude difference between the station and the typhoon center, , are the latitudes of the tide station and the typhoon center, respectively. Further, the segmented constant is used to smooth the typhoon grade data, the missing time is filled with the previous original value, and for the near-center maximum wind speed and the central pressure, the dynamic window Savitzky-Golay filter and the segmented linear interpolation are used for time series preprocessing, and the window length of the filter is set to wherein is a half window length, and is dynamically taken as 1-3, is an odd number of 3-7, and the order of the fitting polynomial is 2, and the filter formula is as follows: ， Wherein, t is the current time, i is the discrete step number of the current time t, is a filter coefficient, is the original data in the window, is the filtered data; For the filtered 1-hour data, piecewise linear interpolation is used to generate 15-minute data, and the time data are set as , , , , , The interpolation result at time is: For the residual series of water level, cubic spline interpolation is used to improve the resolution, and the interval is set as , m is the number of adjacent intervals in the series, and the cubic polynomial is constructed as , , , , , ; the coefficients are obtained by solving the linear equations , , , , and the 15-minute resolution residual series of water level is generated.

4. The regional storm surge forecasting method of claim 1, wherein: The specific steps for constructing a dynamic edge weight graph structure based on a spatiotemporal causal correlation analysis method are as follows: a dynamic weight graph structure is constructed based on Liang-Kleeman information flow, first, define a graph G=(V,E,W), wherein V={A,B,C,D,E} is a tide gauge node, E is a causal edge between nodes, and W is a weight matrix of the edge; secondly, the information flow intensity between stations is calculated, the Liang-Kleeman information flow theory is used to quantify the causal influence on the storm surge residual time series of each station, and the formula is as follows: , where A is a driving variable in the time series, B is a response variable in the time series, is a time derivative of A, is a covariance, is a variance; Introducing a time lag parameter τ, we calculate the information flow intensity under different τ values ​​and select one that makes the information flow intensity more stable. The largest τ is used as the optimal lag, and causal relationships are filtered through information flow values.

5. The regional storm surge forecasting method of claim 4, wherein: The specific steps for constructing a dynamic edge weight graph structure based on a spatiotemporal causal correlation analysis method also include: a graph neural network based on GATv2 is constructed, the features of the node s and the adjacent nodes g are input, the attention coefficient is calculated through linear transformation and LeakyReLU activation function, and the formula is as follows: ， ， wherein, is the unnormalized attention score, is the attention vector, T is the matrix transpose symbol, is a trainable weight matrix, , is the feature vector of node s, g, is the neighborhood of node s, is the attention coefficient, softmax is the normalization function, exp is the exponential function, and l denotes the loop index.

6. The regional storm surge forecasting method of claim 1, wherein: The specific steps for fusing the graph attention mechanism of GATv2 and the temporal feature extraction capability of GRU are as follows: the features of the adjacent nodes are weighted and aggregated by the attention coefficient, and an activation function is applied, the spatial representation of the node is generated by averaging the outputs of the multiple attention heads, and a GRU sequence model is constructed; the reset gate and the update gate are defined, the gate value and the candidate state are obtained according to the current input and the hidden state at the last moment, the state update is completed by element gating, the hidden state is calculated, the hidden state of the last layer of GRU is taken as input to a linear layer to map to a prediction value, and the storm surge prediction value is output for each station, and the formula is as follows: ， wherein is an output layer weight matrix, is an output layer bias term, is a hidden state vector at time step T, is a prediction.

7. The regional storm surge forecasting method of claim 6, wherein: The specific steps of fusing the graph attention mechanism of GATv2 with the time sequence feature extraction capability of GRU include: The parameters are set for model training, the mean square error is used as the loss function, the Adam optimizer is used to update the model parameters, the model overfitting is evaluated by the validation set to adjust the number of GATv2 attention heads and the number of GRU layers, and the trained model is applied to the independent test set to generate the preliminary storm surge prediction results.

8. The regional storm surge forecasting method of claim 1, wherein: The specific steps of model evaluation and result optimization are: The root mean square error, mean absolute error and determination coefficient are calculated, the consistency of the predicted sequence and the measured sequence is evaluated, the differences in the above indicators of TGAT model, SVR, LSTM and atLSTM are compared under the same data and evaluation process, the horizontal comparison results of different models are obtained, the spatial and temporal distribution characteristics of the prediction residual are analyzed by statistics and visualization, the high error interval is located, and the input data quality is checked in parallel. The time neighborhood of the sample determined as the high error sample is modified by using the local regression smoothing technique.

Images