A nearshore wave intelligent forecasting method, medium and system suitable for digital simulation

Abstract

The application provides a nearshore wave intelligent prediction method, medium and system suitable for digital simulation, and belongs to the technical field of marine prediction. The application constructs a typhoon event list by collecting historical typhoon paths, carries out large-scale simulation by driving a numerical wave mode with reanalysis wind field, drives a nearshore high-resolution non-hydrostatic mode with the output two-dimensional sea wave spectrum as a boundary condition, constructs a nearshore high-precision sea wave data set, screens time sample points with effective wave height exceeding a threshold value from the data set, establishes an end-to-end mapping data set with a boundary water level sequence as input and a full-field wave surface height program sequence as output, continuously optimizes the model through hierarchical precision evaluation and high wave height sample supplement mechanism, realizes rapid prediction of the full-field wave surface height time sequence of nearshore typhoon waves, and solves the technical problems that the cost of high-precision numerical simulation of nearshore typhoon waves is too large and the near-real-time prediction time limit of a digital simulation system cannot be met under extreme sea conditions.

Description

Technical Field

[0001] This invention belongs to the field of marine forecasting technology, and more specifically, relates to a method, medium, and system for intelligent forecasting of nearshore waves suitable for digital simulation. Background Technology

[0002] Nearshore typhoon wave forecasting is a core requirement for marine disaster prevention and mitigation. Traditional methods primarily rely on numerical models such as SWAN and SWASH, which simulate nearshore wave propagation, breaking, and flooding processes driven by typhoons with high precision by solving wave action conservation equations and shallow-water non-hydrostatic equations. However, when these numerical models have a grid resolution on the order of 10 meters and cover a computational area of ​​several square kilometers near the coast, a single simulation often takes several hours or even longer, making it difficult to meet the near-real-time forecasting requirements of digital simulation systems. In digital twin coastal protection scenarios, when a typhoon's path deviates rapidly or its intensity increases sharply, the system needs to complete a full-field wavefront elevation time-series forecast within a minute-level time window. Traditional numerical models, due to the iterative convergence mechanism of solving partial differential equations, exhibit a non-linear increase in computational cost with increasing spatial resolution, particularly under high-wave and extreme sea conditions. In other words, existing technologies suffer from the technical problem of excessive computational cost for high-precision numerical simulation of nearshore typhoon waves and the inability to meet the near-real-time forecasting requirements of digital simulation systems under extreme sea conditions. Summary of the Invention

[0003] In view of this, the present invention provides a method, medium and system for intelligent forecasting of nearshore waves suitable for digital simulation, which can solve the technical problems in the prior art where the numerical simulation calculation of nearshore typhoon waves takes too long and is difficult to meet the near real-time forecasting timeliness requirements of digital simulation systems.

[0004] The present invention is implemented as follows: The first aspect of the present invention provides a method for intelligent nearshore wave forecasting suitable for digital simulation, comprising the following steps:

[0005] Collect historical typhoon track, intensity, and timing data, filter typhoon events covering different track types and intensity levels, and construct a list of typhoon events;

[0006] Driven by reanalysis wind field, based on global topographic and water depth data, a nearshore wave simulation model is used to conduct large-scale wave numerical simulation of the selected typhoon process, and output the significant wave height, mean period, peak period, mean wave direction and two-dimensional wave spectrum of the target point.

[0007] Using the two-dimensional wave spectrum output by the nearshore wave simulation model as the boundary condition, a high-resolution and refined simulation of the nearshore area is carried out using a shallow water non-hydrostatic numerical model. The simulation outputs wave surface elevation, current velocity, significant wave height, wave rise and average current velocity, and a high-precision nearshore wave dataset is constructed.

[0008] From the nearshore high-precision ocean wave dataset, samples of moments when the effective wave height exceeds the screening threshold are selected, and the training set, validation set and test set are divided proportionally. Using the boundary water level time series as input and the full-field wave surface elevation time series as output, an end-to-end mapping dataset is constructed.

[0009] Based on the end-to-end mapping dataset, a wavefront time series prediction model is trained. The wavefront time series prediction model consists of a U-ConvGRU branch coupled with a U-Net branch. The U-ConvGRU branch is responsible for wavefront sequence prediction at large-scale time steps, while the U-Net branch is responsible for wavefront sequence reconstruction at small-scale time steps. The two branches are connected in series to output a complete wavefront elevation time series.

[0010] The trained wavefront time series forecast model is used to infer the test set. The effective wave height is calculated from the wavefront height sequence. The forecast accuracy is evaluated according to the effective wave height. When the correlation coefficient of the high wave height interval is lower than the correlation coefficient threshold, weighted oversampling is applied to the high wave height samples and residual connection and attention mechanism are introduced for retraining.

[0011] The steps for constructing a list of typhoon events specifically involve filtering typhoon events within the calculated area based on typhoon path data published in the typhoon path database, covering different landfall locations, durations of impact, and intensity levels, with a minimum of 20 typhoons selected.

[0012] The reanalysis wind field is the fifth-generation global atmospheric reanalysis product of the European Centre for Medium-Range Weather Forecasts, with a horizontal resolution of 0.25° and a temporal resolution of 1 hour. The global topographic and water depth data has a resolution of 1 arc minute. The nearshore wave simulation model uses an unstructured grid with a time step of 10 minutes and an output step of 1 hour.

[0013] The governing equations of the nearshore wave simulation model are wave action conservation equations with source and sink terms. The source and sink terms consist of wind intake, white wave dissipation, bottom friction dissipation, and wave-wave nonlinear interaction terms. Bottom friction dissipation adopts the eddy viscosity model, the fourth-order wave-wave interaction adopts the discrete spectrum interaction approximation method, and the third-order wave-wave interaction adopts the ensemble third-order approximation method.

[0014] The governing equations of the shallow water non-hydrostatic numerical model include a continuity equation and a momentum equation. The spatial step size is 10m in both the x and y directions. The boundary condition type is a weak reflection boundary. The bottom friction adopts the Manning coefficient. The turbulence closure adopts the Smaglinsky model. The wave breaking adopts the still water front approximation method. The breaking parameters are determined by the breaking threshold parameter and the breaking recovery parameter.

[0015] Specifically, the step of constructing the end-to-end mapping dataset involves taking the first row of the water surface elevation data output by the shallow water non-hydrostatic numerical model as the boundary water level time series, and the subsequent rows as the full-field water level distribution. The screening threshold is 1.0m, the ratio of the training set, validation set and test set is 8:1:1, and the data normalization adopts the channel-wise maximum-minimum normalization processing.

[0016] The input of the U-ConvGRU branch includes a boundary water level sequence and a two-dimensional terrain data tensor. The one-dimensional boundary water level sequence is padded and expanded into a two-dimensional tensor consistent with the terrain data. After being fused with the terrain data in the channel dimension, it is fed into the encoder. The encoder consists of multiple convolutional gated recurrent unit blocks, and the decoder consists of multiple skip convolutional gated recurrent unit blocks. The hidden state of the corresponding layer of the encoder is introduced into the decoder through skip connections.

[0017] The U-ConvGRU branch introduces a multi-scale wavelet transform layer to decompose the input boundary water level sequence into high-frequency detail signals and low-frequency trend signals, which are then sent to different channels. A gated routing weight mechanism is used at the neuron connection to dynamically trigger the hidden layer backtracking based on the real-time prediction error rate of the streaming data, and to perform local weight correction on time steps where the accumulated error exceeds the error threshold.

[0018] Among them, the internal weight update of the wavefront time series prediction model introduces an ant colony pheromone enhancement mechanism, which regards the absolute value of the weight gradient between neurons in each layer as the pheromone concentration. The weights of high contribution paths are increased proportionally after each iteration, while the weights of low contribution paths are decayed according to the evaporation rate. When the pheromone concentration is lower than the evaporation threshold, the corresponding weights are blocked from participating in the update in the next iteration.

[0019] Among them, the memory allocation adopts the Gherkin graph dynamic allocation algorithm, which uses the feature map size and channel number of each convolutional layer as the graph node weights and the inter-layer data flow as the directed edges. The shortest path algorithm is used to solve the storage and release order of inter-layer activation values ​​with the minimum peak memory usage. The memory allocation adopts the time step circular buffer mechanism, which only retains the boundary water level data of the current time step and the previous few time steps.

[0020] The U-Net branch adopts a four-layer encoder-decoder structure, introduces skip connections to fuse features at different scales, uses group normalization instead of batch normalization in the normalization layer, and uses a linear rectified function as the activation function. It is responsible for reconstructing the small-scale wavefront sequence between adjacent large-scale prediction times of the U-ConvGRU branch.

[0021] The inter-layer data loop adopts a streaming pipeline parallel approach, and the encoder downsampling layer, decoder upsampling layer and output branch form a three-stage pipeline in the time step dimension; the wavefront time series prediction model training adopts distributed data parallel training, the loss function adopts mean square error, and the optimizer adopts the adaptive moment estimation optimization algorithm.

[0022] The effective wave height is calculated from the wavefront height sequence using Fast Fourier Transform, the one-sided spectrum is estimated by the discrete Fourier transform coefficients of the wavefront discrete sequence and the sampling time interval, and the effective wave height is obtained by spectral integration. The correlation coefficient threshold is determined by statistically analyzing the correlation coefficient on the validation set when it is lower than the correlation coefficient benchmark value for several consecutive rounds.

[0023] Among them, the wave height classification dynamic adjustment function adjusts the large-scale prediction step interval of the U-ConvGRU branch based on the current full-field average effective wave height. The full-field average effective wave height is calculated in real time by the full-field wavefront height program sequence of the previous output time through fast Fourier transform and is then input into the step size control module of the U-ConvGRU branch.

[0024] The screening threshold is 1.0m; the correlation coefficient baseline value is 0.70, and the number of consecutive rounds is 5; the evaporation rate ranges from 0.05 to 0.3; the large-scale prediction step interval is 10s when the average significant wave height is less than 1.0m, 7s in the range of 1.0m to 2.0m, 5s in the range of 2.0m to 3.0m, and 3s when it is not less than 3.0m; the search range of the weighted oversampling coefficient is an integer between 1 and 10.

[0025] A second aspect of the present invention provides a computer-readable storage medium storing program instructions that, when executed in a computer, perform the aforementioned intelligent nearshore wave forecasting method suitable for digital simulation.

[0026] A third aspect of the present invention provides a nearshore wave intelligent forecasting system suitable for digital simulation, comprising the aforementioned computer-readable storage medium, wherein the system is a computer, the computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor for executing program instructions stored in the computer-readable storage medium.

[0027] This invention constructs a wavefront time-series forecasting model based on the coupling of U-ConvGRU and U-Net, using high-precision numerical simulation results of historical typhoon processes as training data. It establishes an end-to-end mapping relationship from boundary water level time series to full-field wavefront elevation time series, solving the technical problem that traditional numerical simulations are too time-consuming under extreme sea conditions and cannot meet the near-real-time forecasting requirements. The U-ConvGRU branch introduces multi-scale wavelet transform to decouple high-frequency breaking signals from low-frequency trend signals, enabling the model to capture long-period wave evolution patterns while preserving near-shore high-frequency nonlinear details. A gated routing weight mechanism suppresses error accumulation in long-term forecasts through real-time error feedback. An ant colony pheromone enhancement mechanism continuously strengthens high-contribution weight paths, ensuring the model's forecasting accuracy in high-wave, low-sample intervals. The synergistic effect of the Gherkin diagram memory allocation algorithm and the streaming pipeline parallel mechanism ensures that the inference speed meets the near-real-time requirements of digital simulation systems. In summary, this invention solves the technical problems mentioned in the background art, such as the excessive computational cost of high-precision numerical simulation of nearshore typhoon waves, which makes it impossible to meet the near-real-time forecasting timeliness requirements of digital simulation systems under extreme sea conditions. Attached Figure Description

[0028] Figure 1 This is a flowchart of the method of the present invention.

[0029] Figure 2 A comparison chart of the accuracy error indices for forecasting different effective wave heights.

[0030] Figure 3 This is a comparison chart of the boundary water level sequence during the peak period of the typhoon and the time series forecast of the wave surface elevation of the entire field. Detailed Implementation

[0031] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below.

[0032] like Figure 1 The diagram shown is a flowchart of a nearshore wave intelligent forecasting method suitable for digital simulation provided by the first aspect of the present invention. This method includes the following steps:

[0033] S01. Collect data on the paths, intensities, and times of typhoons passing through the target area of ​​the Bohai Sea from 1975 to 2025, select no fewer than 20 typhoon events covering different path types and intensity levels, and construct a list of typhoon events.

[0034] S02. Driven by the ERA5 reanalysis wind field, based on ETOPO1 global topographic and water depth data, the SWAN model is used to carry out large-scale wave numerical simulation of the typhoon screening process, and output the significant wave height, average period, peak period, average wave direction and two-dimensional wave spectrum of the target point.

[0035] S03. Using the two-dimensional wave spectrum output by the SWAN mode as the boundary condition, the SWASH mode is used to carry out high-resolution and refined simulation of the nearshore area, outputting five types of elements: wave surface elevation, current velocity, significant wave height, wave rise and average current velocity, and constructing a high-precision nearshore wave dataset.

[0036] S04. Select moments with effective wave height greater than 1.0m from the nearshore high-precision ocean wave dataset, and obtain a total of 484 moment samples. Divide the dataset into training set, validation set and test set in a ratio of 8:1:1. Use the boundary water level time series as input and the full field wave surface elevation time series as output to construct an end-to-end mapping dataset.

[0037] S05. Based on the end-to-end mapping dataset, train the wavefront time series prediction model. The model is composed of a U-ConvGRU branch and a U-Net branch coupled together. The U-ConvGRU branch is responsible for wavefront sequence prediction at large-scale time steps, and the U-Net branch is responsible for wavefront sequence reconstruction at small-scale time steps. The two branches are connected in series to output a complete wavefront elevation time series of 2400 time steps.

[0038] S06. The trained wavefront time series forecast model is used to infer the test set. The effective wave height is calculated from the wavefront height sequence using fast Fourier transform. The mean absolute error, mean relative error, root mean square error, and correlation coefficient are statistically analyzed. The forecast accuracy is evaluated according to the effective wave height greater than 1m, greater than 2m, and greater than 3m. When the correlation coefficient in the effective wave height greater than 3m is lower than the threshold, the high wave height sample supplementation and model structure optimization process is triggered. Weighted oversampling is applied to the high wave height samples, and residual connections and attention mechanisms are introduced for retraining.

[0039] The typhoon event list was selected based on typhoon path data published by the China Typhoon Network. The typhoon paths were required to pass through the calculation area of ​​117° to 127°E and 31° to 41°N, covering different landfall locations, duration of impact, and intensity levels, in order to ensure the dataset's ability to represent various sea conditions. A total of 22 typhoons were selected between 1975 and 2025.

[0040] The ERA5 reanalysis data is the fifth-generation global atmospheric reanalysis product of the European Centre for Medium-Range Weather Forecasts, with a horizontal resolution of 0.25° and a temporal resolution of 1 hour, covering the period from 1959 to the present. It is provided in NetCDF format and includes hourly atmospheric, land and ocean climate variables. The wind field data used in this study to drive the SWAN model covers the period from 2002 to 2021.

[0041] Among them, ETOPO1 is the highest resolution global topographic relief data known to date, which includes two types of data: land topography and ocean depth. It has a resolution of 1 arcminute and has two versions: ice surface and bedrock. This study uses its ocean depth data to provide water depth topographic conditions for the SWAN model.

[0042] The SWAN model, or nearshore wave simulation model, was developed by Delft University of Technology in the Netherlands based on the third-generation wave model WAM. It is suitable for simulating the propagation of random short-peaked waves in nearshore waters. The SWAN model uses an unstructured grid with a grid range of 117° to 127°E and 31° to 41°N, a time step of 10 minutes, and an output step of 1 hour. The physical processes include linear growth, exponential growth, white crown dissipation, four-wave nonlinear interaction, three-wave nonlinear interaction, and bottom friction. The governing equations of the SWAN model are wave action conservation equations with source and sink terms, in Cartesian coordinate form. ; spherical coordinate system form is ;The amount of wave action in the formula With energy spectral density satisfy , For phase frequency, For wave direction, and For the speed of propagation in geospatial space, For frequency space propagation speed, The speed at which the wave propagates into space. Energy source / sink term. Wind intake White Wave Dissipation Item Bottom friction dissipation term Nonlinear interaction term with wave Composition, that is The wind intake term is calculated using the sum of linear and exponential growth: In the formula It is a linearly growing term. The term is exponentially increasing, determined by the wave's frequency, direction, wind speed, and wind direction. White wave dissipation is modeled using the Hasselmann (1974) pulse principle. In the formula For wave number, and These are the average frequency and the average wavenumber, respectively. The coefficient is related to wave steepness. Bottom friction dissipation is calculated using: In the formula When the coefficient of friction is low, and the eddy viscosity model of Madsen et al. (1988) is used, , The velocity of the water particle trajectory. The bottom friction factor is taken as the equivalent roughness length of the seabed. hour The values ​​are approximately 0.01 to 0.02. Fourth-order wave-wave interactions were calculated using the discrete spectrum interaction approximation method of Hasselmann et al. (1985), and third-order wave-wave interactions were calculated using the ensemble third-order approximation method of Eldeberky et al. (1996). SWAN model validation was performed using the process of Typhoon Lekima (1909), with comparative validation based on measured significant wave height and mean period data from the Xiaomaidao and Qianliyan stations.

[0043] The SWASH model is a shallow-water non-hydrostatic numerical model suitable for nearshore to coastal wave simulation, and it has good applicability to wave rise and wave-generated currents caused by wave breaking under complex terrain. The governing equations of the SWASH model include continuity equations and momentum equations: In the formula For time, This refers to the water surface elevation relative to the still water surface. and They are respectively and Vertical average velocity For the total water depth, Still water is deep. It is the acceleration due to gravity. Non-hydrostatic pressure, The coefficient of friction is based on the Manning roughness coefficient. , , , This represents the horizontal turbulent stress term. In the free surface boundary condition, the vertical velocity at the free surface is... satisfy Bottom vertical flow velocity Local mass conservation is The free surface pressure boundary condition is: The bottom horizontal velocity gradient satisfies The solid wall adopts a total internal reflection boundary condition. The outflow adopts a radiation boundary condition. In the formula For wavefront, For wave speed, take SWASH mode space step size and All directions were 10m, the total computation time was 1200s, the initial time step was 0.1s, a dynamic time step was used, the boundary condition type was weak reflection boundary, the bottom friction used the Manning coefficient, the turbulence closure used the Smagorinsky (1963) model, and the wave breaking used the still water front approximation method of Smit et al. (2013). Breaking parameters... , The porosity was set to 1 for the marine region and 0.75 for the land region. The output statistics duration was set to 600 seconds after the simulation. The number of boundary input wave spectra was two, located at 450m and 2200m from the origin, respectively. The output elements included five categories: water surface elevation, current velocity, significant wave height, wave rise, and average current velocity.

[0044] When constructing the end-to-end mapping dataset, the SWASH mode outputs water surface elevation data with the following dimensions: The simulation duration is 1200s, and the output time step is 0.5s. ; Retrieve the first row of data As a time series of boundary water levels, the last 93 rows of data As a whole field water level distribution; a total of 484 times with effective wave height greater than 1.0m were screened. After removing samples with a length of less than 2400 time steps, the training set consisted of 378 samples, the validation set consisted of 48 samples, and the test set consisted of 46 samples.

[0045] The wavefront time-series prediction model is composed of a coupled U-ConvGRU branch and a U-Net branch. The specific structure of the wavefront time-series prediction model is as follows: the model adopts an encoder-decoder topology, with a U-ConvGRU large-scale time branch and a U-Net small-scale time branch coupled in series. The input to the U-ConvGRU branch contains two types of data: one is the time window. to The boundary water level sequence, with dimension indivual One-dimensional time series array, two-dimensional The two-dimensional terrain data tensor, after being padded and expanded into a two-dimensional tensor consistent with the terrain data by the one-dimensional boundary water level sequence, is then fused with the terrain data in the channel dimension. The fused data has the following dimensions: The input is fed into the encoder. The encoder consists of multiple ConvGRU blocks, with a basic structure of convolutional layers, Leaky ReLU activation functions, and ConvGRU layers. In the second and third blocks, average pooling layers are added after the activation functions for downsampling, gradually extracting features and reducing spatial dimensionality. The decoder consists of multiple Skip-ConvGRU blocks, with a basic structure of Skip-ConvGRU layers, deconvolutional layers, and Leaky ReLU activation function layers. Skip connections introduce the hidden states of the corresponding layers of the encoder into the decoder, gradually restoring spatial resolution. The topmost block replaces the deconvolutional layer with a convolutional layer. After decoding, the output branch is connected, consisting of multiple convolutional layers and SiLU activation functions stacked together, post-processing the decoded features into the spatial distribution of water levels in the target area. A multi-scale wavelet transform layer is introduced into the model to decompose the input boundary water level sequence into high-frequency detail signals and low-frequency trend signals, which are fed into different channels, enabling the encoder to simultaneously capture the instantaneous fragmentation details and long-term trend evolution of wave propagation. A gated routing weight mechanism is employed at neuron connections, dynamically triggering hidden layer backtracking based on the real-time prediction error rate of streaming data. Local weight correction is performed at time steps where the accumulated error exceeds a threshold. The decoder is designed with a multi-level residual feedback structure, which enhances the ability to recover local wave-breaking details by nonlinearly concatenating low-level fine-grained spatial features with high-level abstract temporal features. The internal weight update of the model introduces an ant colony pheromone enhancement mechanism, treating the weight gradient between neurons in each layer as pheromone concentration. The weights of high-contribution paths are increased proportionally after each iteration, while the weights of low-contribution paths are increased according to the pheromone evaporation rate. (Values ​​range from 0.05 to 0.3, determined through multiple rounds of iterative experiments on the validation set, with the minimum validation loss as the criterion) Attenuation is applied to dynamically remove failed feature maps. The memory allocation strategy is based on the Gehrler graph dynamic allocation algorithm, using the feature map size and number of channels of each convolutional layer as graph node weights. Using the inter-layer data flow as directed edges, the shortest path corresponding to the memory allocation scheme is calculated, ensuring that the memory usage of each layer's feature map meets the constraint that the peak value does not exceed the single-card capacity. Memory allocation employs a time-step circular buffer mechanism, retaining only the current time step and previous time steps. Boundary water level data for each time step is used, and the remaining historical data is released immediately after inference to reduce peak memory usage. Inter-layer data looping adopts a streaming pipeline parallel approach, with the encoder downsampling layer, decoder upsampling layer, and output branch forming a three-stage pipeline in the time step dimension to reduce inter-layer waiting latency. The U-Net branch adopts a four-layer encoder-decoder structure, introducing skip connections to fuse features at different scales. The normalization layer uses GroupNorm instead of BatchNorm, and the activation function is ReLU, which is responsible for reconstructing the small-scale wavefront sequence between adjacent large-scale prediction times of the U-ConvGRU branch. The number of output time steps within each large-scale interval is the large-scale step size divided by 0.5s minus 1 time step. The technical effects of the wavefront time-series forecasting model are as follows: the multi-scale wavelet transform layer decouples and processes high-frequency breaking signals and low-frequency trend signals separately, enabling the encoder to capture the long-period wave evolution law without losing high-frequency nonlinear details in nearshore shallow waters; the gated routing weight mechanism triggers hidden layer backtracking through real-time error feedback, effectively suppressing the decline in forecast accuracy caused by the gradual accumulation of errors in long-term time-series forecasts; the ant colony pheromone enhancement mechanism continuously strengthens high-contribution weight paths during iteration, while low-contribution paths automatically decay and exit, allowing the model to maintain effective capture of key wave high evolution paths even in the limited sample range of extreme sea states during typhoons; the synergistic effect of the Göhrl diagram memory allocation algorithm and the streaming pipeline parallel mechanism ensures balanced resource usage of feature maps and data loops in each layer during multi-card distributed training, avoiding training interruption caused by single-card memory overflow, and significantly improving the overall inference speed compared to serial processing, meeting the timeliness requirements of digital twin systems for near-real-time forecast updates. The steps for establishing the training dataset for the wavefront time series prediction model specifically include: based on the wavefront elevation time series of 2400 time steps output by SWASH mode, extracting the first row of boundary water level sequences as input samples, and the following 93 rows of full-field wavefront elevations as corresponding labels; filtering 484 time-step samples with effective wave heights greater than 1.0m; removing data with a time series length of less than 2400 steps; and dividing the dataset into a training set of 378 sets, a validation set of 48 sets, and a test set of 46 sets in an 8:1:1 ratio; and using channel-wise maximum-minimum normalization for data normalization. The steps for training the wavefront time series prediction model specifically include: using an Intel Xeon Platinum 8488C processor and an NVIDIA L20×4 graphics card as the training platform, employing distributed data parallel training; using mean squared error as the loss function; using Adam as the optimizer; and initializing the learning rate through learning rate range testing experiments on the training set (from...). to The interval scanning is used to determine the value corresponding to the fastest rate of loss decrease in the validation set. The U-ConvGRU branch and the U-Net branch are first trained independently until the validation loss converges, and then jointly fine-tuned. Weighted oversampling is applied to high-wavelength samples with effective wavelength greater than 3m. The weight coefficients are determined by performing a grid search on the validation set with the goal of maximizing the correlation coefficient of the high-wavelength interval. The search range is an integer between 1 and 10.

[0046] ConvGRU, or Convolutional Gated Recurrent Unit, replaces the fully connected operations in GRU with convolutional operations. While retaining the gating mechanism's ability to selectively remember and forget temporal information, it captures spatial features by leveraging the local receptive field characteristics of convolution. The core components include update gates, reset gates, and candidate memory states, and all component operations are implemented through convolutional operations.

[0047] Skip-ConvGRU is an improved unit that introduces skip connections on the basis of ConvGRU. It introduces the hidden state of the corresponding layer of the encoder into the timing unit of the decoder in a skip connection manner, preserving the fine-grained spatial features captured by the encoder and enhancing the decoder's ability to reconstruct local details.

[0048] Among them, the Gehrler graph dynamic allocation algorithm is a memory scheduling method based on graph theory. It uses the product of the feature map size and the number of channels of each layer of the neural network as the weight of the graph nodes, and uses the forward propagation and backward gradient flow between layers as directed edges to construct a directed weighted graph. The algorithm uses Dijkstra's shortest path algorithm to solve the storage and release order of inter-layer activation values ​​that minimizes peak memory usage, thereby reducing the peak memory usage without affecting the correctness of gradient calculation.

[0049] In this mechanism, the ant colony pheromone enhancement treats the absolute value of the weight gradient between neurons as the pheromone release amount, and after each backpropagation, the cumulative pheromone concentration of each element in the weight matrix is ​​adjusted accordingly. Update, in the formula For the first Neuron to the first The concentration of pheromones in neuronal connections For evaporation rate, This represents the absolute value of the corresponding weight gradient in this backpropagation. When the evaporation rate falls below a preset evaporation threshold, the corresponding weight is blocked from participating in the update in the next iteration, thereby dynamically eliminating failed feature mapping paths; evaporation rate The value ranges from 0.05 to 0.3. The optimal value is determined by conducting multiple rounds of grid search experiments on the training set, with the maximum correlation coefficient in the interval where the effective wave height is greater than 3m as the criterion.

[0050] The significant wave height is calculated by fast Fourier transform from the wavefront height sequence, and the single-sided spectrum estimation formula is as follows: In the formula , , For sample size, The sampling time interval, For the centered wavefront discrete sequence The Discrete Fourier Transform of the 1st There are complex coefficients, and the covariance function is... Discretization estimate is The effective wave height is obtained by spectral integration, and the spectrum... With covariance function The relationship is .

[0051] Among them, the method for determining the correlation coefficient threshold for the effective wave height greater than 3m in the error assessment grading threshold is as follows: statistically analyze the correlation coefficient change curve of different effective wave height intervals with training rounds on the validation set, and take the value corresponding to the correlation coefficient being lower than 0.70 within 5 consecutive rounds as the threshold for triggering the high wave height sample supplementation and model structure optimization process. The specific value of 0.70 is determined by statistically analyzing the prediction residual distribution of samples with effective wave height greater than 3m in the validation set, with the constraint that the proportion of samples with residuals exceeding twice the standard deviation does not exceed 20%, and through multiple rounds of iterative experiments.

[0052] Among them, the wave height classification dynamic adjustment function Based on the current average effective wave height of the entire field (Unit: m) Calculation of the large-scale prediction step size interval used to adjust the U-ConvGRU branch. (Unit: s): When hour, ;when hour, ;when hour, ;when hour, The grading thresholds of 1.0m, 2.0m, and 3.0m correspond to the statistical intervals for the grading of significant wave height errors. The values ​​for each level were determined through enumeration experiments on the validation set with the goal of minimizing prediction error. In the experiments, the candidate step size intervals (3s, 5s, 7s, 10s) were trained separately, and the root mean square error of the validation set was statistically analyzed. The value corresponding to the lowest error was taken as the final value. At each inference time, the full-field wavefront height sequence from the previous output time is calculated in real time using a fast Fourier transform and then fed into the step size control module of the U-ConvGRU branch to dynamically adjust the large-scale step size interval of subsequent predictions. This allows for improved inference efficiency with a larger step size under low wave and normal sea state conditions, and ensured forecast accuracy with a smaller step size under high wave and extreme sea state conditions.

[0053] The specific implementation of step S01 is as follows: Historical typhoon data for the target sea area of ​​the Bohai Sea from 1975 to 2025 are collected, including typhoon paths, central pressures, maximum wind speeds, and durations of impact. Data is filtered based on typhoon path data published by the China Typhoon Network, requiring typhoon paths to pass through the calculation area between 117° and 127°E and 31° and 41°N, and covering different landfall locations, durations of impact, and intensity levels to ensure the dataset's ability to represent various sea conditions. A total of 22 typhoons were obtained after filtering, forming a typhoon event list, providing a basis for selecting driving fields for subsequent numerical simulations.

[0054] The specific implementation of step S02 is as follows: Using ERA5 reanalysis data from the European Centre for Medium-Range Weather Forecasts (ECMWF) as the driving wind field, with a horizontal resolution of 0.25° and a time resolution of 1 hour; and using ETOPO1 global topographic and water depth data (resolution 1 arcminute) as the water depth and topographic conditions, a large-scale wave numerical simulation of 22 typhoon events is conducted using the SWAN model (Nearshore Wave Simulation Model). The SWAN model's governing equations are wave action conservation equations with source and sink terms. The wave action density is calculated in spherical coordinates across geographic space, frequency space, and wave direction space. The source and sink terms include wind intake, white wave dissipation, bottom friction dissipation, and wave-wave nonlinear interaction terms. The model uses an unstructured grid, covering 117° to 127°E and 31° to 41°N, with a time step of 10 minutes and an output step of 1 hour. White wave dissipation was performed using the Hasselman (1974) pulse principle model, bottom friction dissipation using the Madsen et al. (1988) eddy viscosity model, and the bottom friction factor was taken as the reference value corresponding to the equivalent roughness length of the seabed of 0.005m. Fourth-order wave-wave interaction was performed using the discrete spectral interaction approximation method, and third-order wave-wave interaction was performed using the ensemble third-order approximation method of Elderbeck et al. (1996). The SWAN model was validated using the process of Typhoon Lekima (1909) as the verification field, and the significant wave height and mean period were compared and verified using the measured significant wave height and mean period at Xiaomaidao and Qianliyan stations. The output elements include significant wave height, mean period, spectral peak period, mean wave direction, and a two-dimensional wave spectrum at the target point.

[0055] The specific implementation of step S03 is as follows: using the two-dimensional wave spectrum of the target point output by the SWAN model as the boundary condition, a high-resolution, refined simulation of the nearshore area is carried out using the SWASH model (shallow water non-hydrostatic numerical model). The governing equations of the SWASH model include a continuity equation and a momentum equation. The momentum equation includes a non-hydrostatic pressure integral term and a bottom friction term based on the Manning roughness coefficient. Turbulence closure adopts the Smaglinsky (1963) model, and wave breaking adopts the still water front approximation method of Schmitt et al. (2013). Breaking parameters... Take 0.6, The time step was set to 0.3. The model space step size was 10m in both the x and y directions, the total computation time was 1200s, the initial time step was 0.1s, a dynamic time step was used, the boundary condition type was weak reflection boundary, the porosity was set to 1 for the ocean region and 0.75 for the land region, the output statistics duration was set to 600s after the simulation, and the number of boundary input wave spectra was 2, located at 450m and 2200m from the origin respectively. The output elements included five categories: wave surface elevation, current velocity, significant wave height, wave rise, and average current velocity, constructing a high-precision nearshore wave dataset.

[0056] The specific implementation of step S04 is as follows: the water surface elevation data output by SWASH mode has the following dimensions: The simulation duration is 1200s, and the output time step is 0.5s. Take the first row of data ( As a boundary water level time series, the last 93 rows of data ( The data represents the overall water level distribution. From the nearshore high-precision wave dataset, moments with a significant wave height greater than 1.0m were selected, resulting in 484 time-step samples. Samples shorter than 2400 time steps were removed, and the dataset was divided into a training set of 378 samples, a validation set of 48 samples, and a test set of 46 samples in an 8:1:1 ratio. Data normalization was performed using channel-wise max-min normalization, forming an end-to-end mapping dataset.

[0057] The specific implementation of step S05 is as follows: The wavefront time series prediction model adopts an encoder-decoder topology, with a U-ConvGRU large-scale time branch and a U-Net small-scale time branch coupled in series. The U-ConvGRU branch input includes the boundary water level sequence within the time window and a two-dimensional terrain data tensor. The one-dimensional boundary water level sequence is padded and expanded into a two-dimensional tensor consistent with the terrain data, and then fused with the terrain data in the channel dimension before being fed into the encoder. The encoder consists of multiple convolutionally gated recurrent unit blocks, with a basic structure of convolutional layers, a leaky linear rectified activation function, and a convolutionally gated recurrent unit layer. Some blocks add an average pooling layer after the activation function for downsampling. The decoder consists of multiple skip convolutionally gated recurrent unit blocks. Through skip connections, the hidden states of the corresponding layers of the encoder are introduced into the decoder to gradually restore the spatial resolution. After decoding, the output branch is stacked with multiple convolutional layers and curve linear unit activation functions. A multi-scale wavelet transform layer is introduced into the model to decompose the input boundary water level sequence into high-frequency detail signals and low-frequency trend signals, which are then fed into different channels. The gated routing weight mechanism dynamically triggers hidden layer backtracking based on the real-time prediction error rate, and performs local weight adjustments at time steps where the accumulated error exceeds the error threshold. The ant colony pheromone enhancement mechanism treats the absolute value of the weight gradient as the pheromone release amount, and adjusts it accordingly. Update pheromone concentration and evaporation rate The value ranges from 0.05 to 0.3. When the pheromone concentration is below the evaporation threshold, the corresponding weight is masked in the next iteration. The Gherkin graph dynamic allocation algorithm uses the product of the feature map size and the number of channels of each convolutional layer as the graph node weight. It solves the inter-layer activation value storage and release order with the minimum peak memory usage using the Dijkstra shortest path algorithm. Memory allocation adopts a time-step circular buffer mechanism, which only retains the boundary water level data of the current time step and the previous few time steps. Inter-layer data circulation adopts a streaming pipeline parallel approach, with the encoder downsampling layer, decoder upsampling layer and output branch forming a three-stage pipeline. The U-Net branch adopts a four-layer encoder-decoder structure, introduces skip connections to fuse features of different scales, uses group normalization for the normalization layer, and uses a linear rectified function for the activation function. It is responsible for reconstructing the small-scale wavefront sequence between adjacent large-scale prediction times of the U-ConvGRU branch. The training platform uses an Intel Xeon Platinum 8488C processor and an NVIDIA L20 quad-card configuration. It employs distributed data parallel training, uses mean squared error as the loss function, and uses an adaptive moment estimation optimization algorithm as the optimizer. The U-ConvGRU branch and the U-Net branch are first trained independently until the verification loss converges, and then jointly fine-tuned.

[0058] The specific implementation of step S06 is as follows: The trained wavefront time-series forecasting model is used to infer the values ​​of 46 samples in the test set. The effective wave height is calculated from the wavefront height sequence using a Fast Fourier Transform. The one-sided spectrum estimation formula is calculated using the discrete Fourier transform coefficients of the wavefront discrete sequence and the sampling time interval. The effective wave height is obtained by integrating the spectrum. The mean absolute error, mean relative error, root mean square error, and correlation coefficient are statistically analyzed, and the forecast accuracy is evaluated according to effective wave heights greater than 1m, 2m, and 3m. When the correlation coefficient in the effective wave height greater than 3m interval is below 0.70 for five consecutive rounds, a high-wave height sample supplementation and model structure optimization process is triggered. Weighted oversampling is applied to the high-wave height samples with effective wave heights greater than 3m. The weight coefficients are determined by performing a grid search on the validation set with the objective of maximizing the correlation coefficient in the high-wave height interval. The search range is integers between 1 and 10. Residual connections and attention mechanisms are introduced for retraining. The wave height classification dynamic adjustment function adjusts the large-scale prediction step interval of the U-ConvGRU branch in real time according to the current mean significant wave height. When the mean significant wave height is below 1.0m, the step interval is 10s; in the range of 1.0m to 2.0m, it is 7s; in the range of 2.0m to 3.0m, it is 5s; and when it is not lower than 3.0m, it is 3s. This allows for a larger step size to improve inference efficiency under low wave height and normal sea state conditions, while a smaller step size to ensure forecast accuracy under high wave height and extreme sea state conditions.

[0059] It should be noted that the key technologies of this invention include: First, the dual-branch serial coupling architecture of U-ConvGRU and U-Net decouples wavefront sequence prediction at large-scale time steps from wavefront sequence reconstruction at small-scale time steps. This allows the model to maintain long-term temporal coverage without losing high-frequency details of nearshore wave breaking. Compared with recurrent networks with a single temporal granularity, this architecture effectively balances computational efficiency and temporal resolution. Second, the synergistic effect of the multi-scale wavelet transform layer and the ant colony pheromone enhancement mechanism. The former decouples the high-frequency breaking signal from the low-frequency trend signal and inputs them separately into the encoder, providing a higher-resolution gradient information basis for the ant colony mechanism. The latter dynamically enhances high-contribution weight paths and attenuates low-contribution paths based on the absolute value of the gradient. The synergy of the two enables the model to maintain effective capture of key wave height evolution paths even in small sample intervals with high waves and extreme sea states, overcoming the limitation of insufficient generalization ability of simple data augmentation methods under extreme sea states. Third, the synergistic effect of the Gehrtu memory allocation algorithm and the streaming pipeline parallel mechanism ensures the balance of multi-card distributed training and inference from the hardware resource scheduling level, so that the overall inference speed meets the near real-time efficiency requirements of the digital simulation system. Compared with the serial processing method, it has a significant advantage in terms of resource utilization.

[0060] It should be noted that when the digital simulation system updates the nearshore wave field in real time, the proportion of high-wave-height moments in the overall wave surface elevation time series is extremely low when the typhoon intensity changes abruptly within a short period. This results in a much higher prediction error for the model in the high-wave-height range compared to normal sea states. This technical problem arises because moments with significant wave heights exceeding 3m under extreme typhoon sea states are already low-probability events in historical simulation data. In the sample set after significant wave height filtering, high-wave-height samples typically account for no more than 10% to 15% of the total sample. This leads to an excessively low weighting of the model's contribution to the loss in the high-wave-height range during training. Weight updates are primarily dominated by a large number of medium- and low-wave-height samples, causing the feature mapping path in the high-wave-height range to gradually degrade due to gradient sparsity. The usual solution to this technical problem is to expand the number of numerical simulation sessions to increase the number of high-wave-height samples. However, due to the time consumption of a single simulation in a near-shore high-resolution non-hydrostatic model, large-scale expansion of simulation sessions is extremely costly in engineering practice. Furthermore, the number of extremely strong typhoons in historical typhoon records is limited, and simply increasing the number of simulation sessions cannot fundamentally change the sparse distribution of high-wave, high-state samples. This invention effectively solves this technical problem. The ant colony pheromone enhancement mechanism treats the absolute value of the weight gradient between neurons as the pheromone release amount. By controlling the automatic decay of low-contribution paths through the evaporation rate parameter, high-contribution weight paths are continuously strengthened during iteration. Even with a limited number of high-wave, high-state samples, the corresponding feature mapping paths can be continuously strengthened during training without being overwhelmed by the gradients of medium- and low-wave, high-state samples. Simultaneously, weighted oversampling is applied to samples with effective wave heights greater than 3m to compensate for the loss function bias caused by uneven sample distribution at the data level. The synergistic effect of these two mechanisms allows the model to maintain high forecast accuracy even in the limited sample range of high-wave, high-extreme sea states without the need for large-scale expansion of numerical simulation sessions.

[0061] A second aspect of the present invention provides a computer-readable storage medium storing program instructions that, when executed in a computer, perform the aforementioned intelligent nearshore wave forecasting method suitable for digital simulation.

[0062] A third aspect of the present invention provides a nearshore wave intelligent forecasting system suitable for digital simulation, comprising the aforementioned computer-readable storage medium. The system is any one of a computer, a server, or a microcontroller. The computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor that executes the program instructions stored in the computer-readable storage medium.

[0063] Specifically, the principle of this invention is:

[0064] The present invention can solve the above-mentioned technical problems because it fundamentally changes the calculation paradigm of nearshore typhoon wave forecasting: traditional numerical models need to resolve the partial differential equation system at each forecast time, and the amount of computation is proportional to the number of grid nodes; while the present invention transforms the results of offline high-precision numerical simulation into training data for neural networks, so that only one forward propagation needs to be performed in the online inference stage, and the amount of computation is related to the number of model parameters after training, decoupled from the grid size, thereby greatly reducing the computation time in the inference stage.

[0065] The technical feasibility of this paradigm shift relies on the following logical chains. First, the cascaded simulation of SWAN and SWASH modes ensures that the training data covers the entire physical process from ocean-scale wave propagation to nearshore breakup and flooding, enabling the neural network to implicitly learn the nonlinear mapping relationship of wave dynamics from the data, rather than relying on manual feature extraction. Second, the encoder-decoder topology of the U-ConvGRU branch, combined with multi-scale wavelet transform, decomposes the boundary water level sequence into high-frequency detail signals and low-frequency trend signals, and inputs them separately. This allows the encoder to simultaneously capture the instantaneous breakup details and long-term trend evolution of wave propagation, avoiding the loss of high-frequency information caused by single-scale feature extraction. Skip connections introduce the hidden states of each layer of the encoder into the corresponding decoder layer, preserving fine-grained spatial features and giving the decoder a strong ability to reconstruct local wave breakup details.

[0066] Secondly, to address the issue of gradual error accumulation in long-term time-series predictions, the gated routing weight mechanism dynamically triggers hidden layer backtracking based on the real-time prediction error rate of streaming data. It performs local weight correction on time steps where the accumulated error exceeds the threshold, thereby blocking the propagation path of error along the time axis from a mechanism perspective.

[0067] Secondly, to address the issue of insufficient model generalization ability in high-wave, high-extreme-sea-state regions due to the scarcity of samples, the ant colony pheromone enhancement mechanism treats the absolute value of the weight gradient as the pheromone release amount. By controlling the automatic decay and exit of low-contribution paths through the evaporation rate parameter, the model can still effectively capture key wave height evolution paths even with limited samples. Simultaneously, weighted oversampling is applied to samples with significant wave heights greater than 3m to compensate for the bias caused by uneven sample distribution at the data level.

[0068] Finally, the Ghertz graph dynamic allocation algorithm optimizes the memory storage and release order of feature maps in each convolutional layer using graph theory's shortest path method. The streaming pipeline parallel mechanism constructs a three-stage pipeline for encoder downsampling, decoder upsampling, and output branching. The synergistic effect of these two mechanisms balances resource usage in multi-GPU distributed training and inference, ensuring, from a hardware perspective, that the inference speed meets the near real-time requirements of digital simulation systems. The wave height-level dynamic adjustment function adaptively adjusts the large-scale prediction step size interval based on the current full-field average significant wave height. Under low wave height and normal sea states, a larger step size improves inference efficiency, while under high wave height and extreme sea states, a smaller step size ensures forecast accuracy, further achieving a balance between accuracy and efficiency.

[0069] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0070] The specific implementation of step S01 is as follows: collect data on the path, intensity, and timing of typhoons passing through the target area of ​​the Bohai Sea from 1975 to 2025; based on the typhoon path data published by the China Typhoon Network, screen typhoon events whose paths pass through the calculation area of ​​117° to 127°E and 31° to 41°N, requiring coverage of different landfall locations, duration of impact, and intensity levels; a total of 22 typhoons were screened, and a list of typhoon events was constructed.

[0071] The specific implementation of step S02 is as follows: Driven by the wind field of the European Centre for Medium-Range Weather Forecasts (ECMWF) fifth-generation global atmospheric reanalysis product ERA5, with a horizontal resolution of 0.25° and a temporal resolution of 1 hour, covering the period from 2002 to 2021. Topographic water depth data uses ETOPO1 global ocean depth data with a resolution of 1 arcminute. A large-scale numerical simulation of waves is conducted on the selected typhoon process based on a nearshore wave simulation model. The governing equation is the wave action conservation equation with source and sink terms, in Cartesian coordinate form:

[0072] ;

[0073] The spherical coordinate system is in the form of:

[0074] ;

[0075] In the formula, Wave action density, in units of , Phase frequency, unit: , Wave direction, unit is , For time, the unit is , and Cartesian coordinates, unit: , and The speed of geographic spatial propagation in the Cartesian coordinate system, in units of , The speed of propagation in frequency space, in units of , The speed at which waves propagate into space, measured in units of 1000 m / s. , Longitude Latitude and These represent the longitudinal and latitudinal propagation velocities in spherical coordinates, respectively, with units of [missing information]. , Energy source sink, unit is All propagation terms on the left side of the equation have the same dimensional structure, while those on the right side... and Dimensions are consistent, all are Wave action With energy spectral density satisfy ,in Units are , Units are , Units are Both sides have the same dimensions. The source and sink terms consist of the following formula:

[0076] ;

[0077] In the formula, For wind intake, For the white wave dissipation term, This is the bottom friction dissipation term. This is the wave-wave nonlinear interaction term, with units of all terms being... The wind intake term is calculated using the sum of linear and exponential growth:

[0078] ;

[0079] In the formula, This is a linearly growing term, with units of... , The exponential growth rate is expressed in units of 1000 ppm. , Energy spectral density, in units of , The product unit is ,and Dimensions consistent and All are determined by the wave frequency, direction, wind speed, and wind direction. White wave dissipation adopts the Hasselman (1974) pulse principle model:

[0080] ;

[0081] In the formula, Wave number, unit: , Average frequency, in units of , The mean wave number is expressed in units of 10 ... , The dimensionless coefficients related to wave steepness, The unit is ,and Dimensions are consistent. The frictional dissipation is:

[0082] ;

[0083] In the formula, When the coefficient of friction is low, using the eddy viscosity model of Madsen et al. (1988):

[0084] ;

[0085] In the formula, The velocity of the water particle trajectory, in units of , It is a dimensionless low friction factor. This is the acceleration due to gravity, in units of 1. ,but Units are , The depth of still water, in units of , It is a dimensionless hyperbolic sine function. Units are , Units are ,and With consistent dimensions, the equivalent roughness length of the seabed is taken. hour The empirical values ​​are approximately 0.01 to 0.02. Fourth-order wave-wave interactions were calculated using the discrete spectral interaction approximation method of Hasselman et al. (1985), and third-order wave-wave interactions were calculated using the ensemble third-order approximation method of Elderbeckie et al. (1996). The model used an unstructured grid with a grid range of 117° to 127°E and 31° to 41°N, a time step of 10 min, and an output step of 1 h. Validation was performed using the process of Typhoon Lekima (1909), comparing the measured significant wave height and mean period data from the Xiaomaidao and Qianliyan stations. Output elements included significant wave height, mean period, spectral peak period, mean wave direction, and a two-dimensional wave spectrum at the target location.

[0086] The specific implementation of step S03 is as follows: using the two-dimensional wave spectrum output by the nearshore wave simulation model as the boundary condition, a high-resolution, refined simulation of the nearshore area is conducted using a shallow water non-hydrostatic numerical model. The governing equations of the shallow water non-hydrostatic numerical model include a continuity equation and a momentum equation. The continuity equation is:

[0087] ;

[0088] The equation for directional momentum is:

[0089] ;

[0090] The equation for directional momentum is:

[0091] ;

[0092] In the formula, This is the elevation of the water surface relative to still water, in units of... , and They are respectively and Vertical average velocity, in units of , Total water depth, in units of , The depth of still water, in units of , This is the acceleration due to gravity, in units of 1. , This is non-hydrostatic pressure, and the unit is... (here) Expressed as kinematic pressure, i.e., pressure divided by the density of water. Vertical coordinates, in units of , The coefficient of friction is based on the Manning roughness coefficient, and its unit is . , Manning's roughness coefficient, in units of , , , , This is the horizontal turbulent stress term, in units of... The units of all terms in the momentum equations are... The dimensions are consistent on both sides. In the free surface boundary condition, the vertical velocity of the free surface is... satisfy:

[0093] ;

[0094] In the formula, Vertical velocity at the free surface, in units of , This refers to the bottom non-hydrostatic pressure (in kinematic pressure form), with units of... , Total water depth, in units of , The vertical flow velocity at the bottom is expressed in units of 1000 m / s. The units of all terms in the equation are The dimensions are consistent on both sides; the vertical flow velocity at the bottom is... satisfy:

[0095] ;

[0096] In the formula, and All units are ,and Dimensionally consistent; local mass conservation is All units are The free surface pressure boundary condition is: The bottom horizontal velocity gradient satisfies The solid wall adopts a total internal reflection boundary condition. ,in The vertical flow velocity at the solid wall is expressed in units of 1 / 2 Å. The outflow adopts a radiation boundary condition:

[0097] ;

[0098] In the formula, Wavefront elevation, in units of , Wave speed, unit: ,Pick The units of all terms in the equation are The dimensions on both sides are consistent. Spatial step size and All directions were 10m, the total calculation time was 1200s, the initial time step was 0.1s, a dynamic time step was used, and the wave breaking was performed using the still water front approximation method of Schmitt et al. (2013). The breaking parameters were... , The porosity was set to 1 for the marine region and 0.75 for the land region. The output statistical duration was set to 600 seconds after the simulation. The number of boundary input wave spectra was 2, located at 450m and 2200m from the origin, respectively. The output elements included five categories: water surface elevation, current velocity, significant wave height, wave rise, and average current velocity, thus constructing a high-precision nearshore wave dataset.

[0099] The specific implementation of step S04 is as follows: The shallow water non-hydrostatic numerical model outputs water surface elevation data with dimensions of... ,in The corresponding simulation duration is 1200s, and the output time step is 0.5s. This represents the spatial grid dimension. Take the first row of data (dimension). As a boundary water level time series, the last 93 rows of data (dimensions) As the overall water level distribution, 484 moments with a valid wave height greater than 1.0m were selected. Samples with a length of less than 2400 time steps were removed. The data were divided into a training set of 378 samples, a validation set of 48 samples, and a test set of 46 samples in an 8:1:1 ratio. Data normalization was performed using channel-wise maximum-minimum normalization to construct an end-to-end mapping dataset.

[0100] The specific implementation of step S05 is as follows: The wavefront time series prediction model is composed of a large-scale time branch (U-ConvGRU branch) and a small-scale time branch (U-Net branch) coupled in series. The convolutional gated recurrent unit replaces the fully connected operations in the gated recurrent unit with convolution operations. The core components include update gate, reset gate, and candidate memory state. The operations of each component are implemented through convolution operations. The update gate... The calculation is as follows:

[0101] ;

[0102] Reset door The calculation is as follows:

[0103] ;

[0104] Candidate memory states The calculation is as follows:

[0105] ;

[0106] Current hidden state The update is as follows:

[0107] ;

[0108] In the formula, It is the sigmoid activation function. For convolution operations, For element-wise multiplication, , , , The convolution kernel weight matrix is... , , For the corresponding bias vector This is the hidden state from the previous moment. Input the feature tensor at the current time step. Hide the current state. and The value range is , The range of values ​​is All of the above quantities are feature tensors, and their dimensions depend on the dimensionless features after input normalization. The input of the U-ConvGRU branch includes a time window. to Boundary water level sequence ( indivual (one-dimensional time series array) and The two-dimensional terrain data tensor, after being padded and expanded into a two-dimensional tensor consistent with the terrain data by the one-dimensional boundary water level sequence, is fused with the terrain data in the channel dimension. The fused data has the following dimensions: Feed into the encoder, This represents the number of channels after fusion. The multi-scale wavelet transform layer will input the boundary water level sequence. Decomposed into high-frequency detail signals Low-frequency trend signals ( (For wavelet decomposition hierarchical indexing), these are fed into different channels, enabling the encoder to simultaneously capture high-frequency nonlinear details and low-frequency trend evolution. The ant colony pheromone enhancement mechanism considers the absolute value of the weight gradient between neurons as the pheromone release amount, and updates the cumulative pheromone concentration of each element in the weight matrix after each backpropagation using the following formula:

[0109] ;

[0110] In the formula, For the first Neuron to the first The concentration of pheromones in neuronal connections The evaporation rate, ranging from 0.05 to 0.3, was determined through multiple rounds of grid search on the validation set, using the highest correlation coefficient in the range of significant wave heights greater than 3 m as the criterion. This represents the absolute value of the corresponding weight gradient in this backpropagation; when When the value falls below a preset evaporation threshold, the corresponding weight is blocked from participating in the update in the next iteration, dynamically eliminating failed feature mapping paths. The Göhrl graph dynamic allocation algorithm uses the product of the feature map size and channel number of each layer of the neural network as the graph node weights, constructing a directed weighted graph with forward propagation and backward gradient flow as directed edges. The Dijkstra shortest path algorithm is used to solve for the storage and release order of inter-layer activation values ​​that minimizes peak memory usage. Wave height graded dynamic adjustment function. Based on the current average effective wave height of the entire field (Unit: m) Adjusting the large-scale prediction step size interval (Unit: s), the rules are as follows: when hour, ;when hour, ;when hour, ;when hour, The values ​​for each level are determined by training on the candidate step size intervals (3s, 5s, 7s, 10s) on the validation set and calculating the root mean square error, and taking the value corresponding to the lowest error. At each inference time step, the high-order sequence of the full-field wavefront from the previous output time step is calculated in real time using a Fast Fourier Transform and then fed into the step size control module of the U-ConvGRU branch. The U-Net branch employs a four-layer encoder-decoder structure, introducing skip connections to fuse features at different scales. The normalization layer uses group normalization instead of batch normalization, and the activation function uses a linear correction unit. This unit is responsible for reconstructing the small-scale wavefront sequence between adjacent large-scale prediction time steps of the U-ConvGRU branch. The number of output time steps within each large-scale interval is... At each moment, the two branches are cascaded to output a complete wavefront elevation time series of 2400 moments. The training platform uses an Intel Xeon Platinum 8488C processor and NVIDIA... The graphics card is used for distributed data parallel training; the loss function is mean squared error, the optimizer uses adaptive moment estimation, and the initial learning rate is obtained by analyzing data from the training set. to The interval scanning is determined; the two branches are first trained independently until the validation loss converges, and then joint fine-tuning is performed; weighted oversampling is applied to high wave height samples with effective wave height greater than 3m, and the weight coefficients are determined by performing a grid search (search range of integers between 1 and 10) on the validation set with the goal of maximizing the correlation coefficient of the high wave height interval.

[0111] The specific implementation of step S06 is as follows: The trained wavefront time series prediction model is used to infer the test set. The effective wave height is calculated from the wavefront height sequence using Fast Fourier Transform. The single-sided spectrum estimation formula is:

[0112] ;

[0113] In the formula, , For the first A discrete angular frequency, in units of (This is the normalized angular frequency, which is dimensionless.) For sample size, The sampling time interval is expressed in units of 1 / 2. , For the centered wavefront discrete sequence The Discrete Fourier Transform of the 1st Each complex coefficient This is a one-sided power spectral density estimate, in units of The right side of the equals sign middle Units are , Units are The product unit is The dimensions on both sides are consistent; the discretized estimate of the covariance function is:

[0114] ;

[0115] In the formula, Time delay The corresponding covariance estimate for each step, in units of , For the first The centered discrete wavefront sequence value at each time step, in units of , The time delay step number; spectrum With covariance function The relationship is:

[0116] ;

[0117] In the formula, Two-sided power spectral density, in units of , It is a continuous covariance function, with units of . , It is a continuous time delay, with units of . , Continuous angular frequency, unit: , It is a dimensionless complex exponential function, and the integral result on the right side is in units of . The dimensions on both sides are consistent, and the significant wave height is obtained by spectral integration. The mean absolute error, mean relative error, root mean square error, and correlation coefficient are statistically analyzed, and the prediction accuracy is evaluated according to significant wave heights greater than 1m, 2m, and 3m. When the correlation coefficient for the significant wave height interval greater than 3m is below the threshold of 0.70 for five consecutive training rounds, a process of supplementing high-wave height samples and optimizing the model structure is triggered. Weighted oversampling is applied to the high-wave height samples, and residual connections and attention mechanisms are introduced for retraining. The threshold of 0.70 is determined through statistical analysis of the prediction residual distribution of samples with significant wave heights greater than 3m in the validation set, with the constraint that the proportion of samples with residuals exceeding twice the standard deviation does not exceed 20%, and through multiple rounds of iterative experiments.

[0118] To better understand and implement this invention, the following is a specific application scenario of this invention, Example 2:

[0119] This embodiment uses a protected section of the Bohai Sea nearshore as a verification scenario. It employs the intelligent nearshore wave forecasting method applicable to digital simulation described in this invention to perform wave surface elevation time-series forecasting for a historical strong typhoon process, and compares and analyzes the results with the corresponding numerical simulation results in the nearshore high-precision ocean wave dataset.

[0120] The typhoon event list was compiled from historical typhoon track data, resulting in 22 typhoons covering different track types and intensity levels across the calculated area between 1975 and 2025. This embodiment selects one typhoon event with a relatively long landfall duration and high intensity as the test sample. During its impact on the target coastline, this typhoon's nearshore significant wave height peaked at over 3.5m, classifying it as a typical high-wave, high-extreme-sea-state event.

[0121] The wind field data, provided at a 1-hour time resolution, was interpolated and used as the driving field for the nearshore wave simulation model. Global topographic and water depth data were nested into the computational domain at a 1-arc minute resolution. The nearshore wave simulation model covered the computational domain with an unstructured grid, outputting two-dimensional wave spectra of significant wave height, mean period, and target point. The model's time step was 10 minutes, and the output step was 1 hour. Using the output two-dimensional wave spectra as boundary conditions, a shallow-water non-hydrostatic numerical model was driven to perform high-resolution simulations of the nearshore area. The spatial step was 10 meters, the output time step was 0.5 seconds, the simulation duration was 1200 seconds, and the output statistical duration was taken from the last 600 seconds of the simulation, resulting in a total of 2400 time steps of wave surface elevation time series.

[0122] From the nearshore high-precision wave dataset, moments with significant wave heights greater than 1.0 m were selected, yielding 484 samples. Among these, 41 samples had significant wave heights greater than 3 m, accounting for approximately 8.5% of all samples, indicating a significant sparsity in the distribution of high-wave-height samples. The training, validation, and test sets were divided in an 8:1:1 ratio, with 378 samples in the training set, 48 in the validation set, and 46 in the test set. Data normalization was performed using channel-wise max-min normalization. The distribution of samples at each wave height level is shown in Table 1.

[0123] Table 1. Distribution of Sample Quantities for Each Significant Wave Height Class

[0124]

[0125] As shown in Table 1, the number of high-wavelength samples with an effective wave height greater than 3m is significantly insufficient. Without weighted oversampling, the loss contribution weight of this interval during model training will be much lower than that of the mid-to-low wavelength interval, leading to a decrease in prediction accuracy for the high-wavelength interval. This invention utilizes the synergistic effect of ant colony pheromone enhancement mechanism and weighted oversampling to continuously strengthen the feature mapping path corresponding to high-wavelength samples during training.

[0126] After the wavefront time-series forecast model was trained, inference was performed on 46 samples from the test set. During the inference phase, the wave height classification dynamic adjustment function adjusted the large-scale prediction step interval of the U-ConvGRU branch in real time based on the overall mean significant wave height. During the peak period of the typhoon process in this test, the overall mean significant wave height exceeded 3.0m, and the step interval automatically switched to 3s, ensuring forecast accuracy under extreme sea conditions. During the recovery period after the typhoon weakened, the overall mean significant wave height dropped below 1.0m, and the step interval switched to 10s, correspondingly improving inference efficiency. Figure 2 As shown, the forecast accuracy evaluation results of each significant wave height classification interval demonstrate the trend of error index changes with significant wave height classification. As the significant wave height classification threshold increases, the root mean square error first increases and then tends to stabilize. This trend is consistent with the physical logic that high wave height samples are compensated by weighted oversampling and the ant colony mechanism strengthens the corresponding feature path.

[0127] The statistical results of each error assessment index are shown in Table 2:

[0128] Table 2. Evaluation Table of Prediction Accuracy for Each Significant Wave Height Class

[0129]

[0130] As shown in Table 2, the correlation coefficient for the interval with significant wave height greater than 3m is 0.91, which is higher than the correlation coefficient threshold of 0.70. This indicates that the retraining process was not triggered, demonstrating that the forecast accuracy of this invention under high wave and extreme sea conditions meets the requirements of a digital simulation system. The average relative error for the interval with significant wave height greater than 3m (6.4%) is lower than that for the interval greater than 1m (8.3%). The physical explanation for this phenomenon is that the high wave and high wave periods correspond to the period of greatest typhoon intensity, a more concentrated spatial distribution of the wave field, and a stable main wave direction, resulting in stronger characteristic regularity captured by the model in this interval. Conversely, the low wave and high wave periods correspond to the complex sea conditions before and after the typhoon's passage, where the superposition of multiple waves increases forecast uncertainty, leading to a slightly higher relative error.

[0131] like Figure 3 As shown in the diagram, the comparison between the boundary water level time series and the corresponding full-field wavefront elevation time series demonstrates the consistency between the full-field wavefront elevation time series output by the wavefront time series forecast model and the output results of the shallow water non-hydrostatic numerical model during the peak period of the typhoon in this test. During the peak period of the typhoon, the instantaneous peak and trough values ​​and phase of the wavefront elevation maintain good agreement with the numerical simulation results, indicating that the multi-scale wavelet transform layer effectively preserves the nonlinear deformation characteristics of nearshore waves in shallow water by decoupling the high-frequency breaking signal and the low-frequency trend signal, rather than only capturing the low-frequency envelope trend.

[0132] While completing the forecast, the computation time of the inference stage of this invention is much lower than that of the simulation time of the shallow water non-hydrostatic numerical model. Under the test hardware environment, the wavefront time series forecasting model performs well for a single sample (2400 time steps). The inference time for a spatial grid simulation is within seconds, while the simulation time for a single shallow-water non-hydrostatic numerical model is on the order of several hours. This difference in computational efficiency stems from the fact that the inference phase only requires forward propagation without iteratively solving partial differential equations. This improvement in computational efficiency allows the digital simulation system to update the entire wavefront elevation time series within a minute-level time window when the typhoon's path or intensity changes, meeting the timeliness requirements of near-real-time forecasts. This solves the technical problem that traditional numerical simulations are computationally expensive under extreme sea conditions and cannot meet the near-real-time forecasting requirements of digital simulation systems.

[0133] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for intelligent nearshore wave forecasting suitable for digital simulation, characterized in that, Includes the following steps: Collect historical typhoon track, intensity, and timing data, filter typhoon events covering different track types and intensity levels, and construct a list of typhoon events; Driven by reanalysis wind field, based on global topographic and water depth data, a nearshore wave simulation model is used to conduct large-scale wave numerical simulation of the selected typhoon process, and output the significant wave height, mean period, peak period, mean wave direction and two-dimensional wave spectrum of the target point. Using the two-dimensional wave spectrum output by the nearshore wave simulation model as the boundary condition, a high-resolution and refined simulation of the nearshore area is carried out using a shallow water non-hydrostatic numerical model. The simulation outputs wave surface elevation, current velocity, significant wave height, wave rise and average current velocity, and a high-precision nearshore wave dataset is constructed. From the nearshore high-precision ocean wave dataset, samples of moments when the effective wave height exceeds the screening threshold are selected, and the training set, validation set and test set are divided proportionally. Using the boundary water level time series as input and the full-field wave surface elevation time series as output, an end-to-end mapping dataset is constructed. Based on the end-to-end mapping dataset, a wavefront time series prediction model is trained. The wavefront time series prediction model consists of a U-ConvGRU branch coupled with a U-Net branch. The U-ConvGRU branch is responsible for wavefront sequence prediction at large-scale time steps, while the U-Net branch is responsible for wavefront sequence reconstruction at small-scale time steps. The two branches are connected in series to output a complete wavefront elevation time series. The trained wavefront time series forecast model is used to infer the test set. The effective wave height is calculated from the wavefront height sequence. The forecast accuracy is evaluated according to the effective wave height. When the correlation coefficient of the high wave height interval is lower than the correlation coefficient threshold, weighted oversampling is applied to the high wave height samples and residual connection and attention mechanism are introduced for retraining.

2. The intelligent nearshore wave forecasting method suitable for digital simulation according to claim 1, characterized in that, The steps to construct a list of typhoon events are as follows: based on the typhoon path data published in the typhoon path database, select typhoon events that pass through the calculation area, covering different landfall locations, duration of impact, and intensity levels, and obtain no fewer than 20 typhoons.

3. The intelligent nearshore wave forecasting method suitable for digital simulation according to claim 2, characterized in that, The reanalysis wind field is the fifth-generation global atmospheric reanalysis product of the European Centre for Medium-Range Weather Forecasts, with a horizontal resolution of 0.25° and a temporal resolution of 1 hour. The global topographic and water depth data has a resolution of 1 arc minute. The nearshore wave simulation model uses an unstructured grid with a time step of 10 minutes and an output step of 1 hour.

4. The intelligent nearshore wave forecasting method suitable for digital simulation according to claim 3, characterized in that, The governing equations of the nearshore wave simulation model are wave action conservation equations with source and sink terms. The source and sink terms consist of wind intake, white wave dissipation, bottom friction dissipation, and wave-wave nonlinear interaction terms. Bottom friction dissipation adopts the eddy viscosity model, the fourth-order wave-wave interaction adopts the discrete spectrum interaction approximation method, and the third-order wave-wave interaction adopts the ensemble third-order approximation method.

5. The intelligent nearshore wave forecasting method suitable for digital simulation according to claim 4, characterized in that, The governing equations of the shallow water non-hydrostatic numerical model include a continuity equation and a momentum equation. The spatial step size is 10m in both the x and y directions. The boundary condition type is a weak reflection boundary. The bottom friction adopts the Manning coefficient. The turbulence closure adopts the Smaglinsky model. The wave breaking adopts the still water front approximation method. The breaking parameters are determined by the breaking threshold parameter and the breaking recovery parameter.

6. The intelligent nearshore wave forecasting method suitable for digital simulation according to claim 5, characterized in that, The steps for constructing the end-to-end mapping dataset are as follows: take the first row of the water surface elevation data output by the shallow water non-hydrostatic numerical model as the boundary water level time series, and take the subsequent rows as the full field water level distribution. The screening threshold is 1.0m, the ratio of the training set, validation set and test set is 8:1:1, and the data normalization adopts the channel-wise maximum-minimum normalization processing.

7. The intelligent nearshore wave forecasting method suitable for digital simulation according to claim 6, characterized in that, The input of the U-ConvGRU branch includes a boundary water level sequence and a two-dimensional terrain data tensor. The one-dimensional boundary water level sequence is padded and expanded into a two-dimensional tensor consistent with the terrain data. After being fused with the terrain data in the channel dimension, it is fed into the encoder. The encoder consists of multiple convolutional gated recurrent unit blocks, and the decoder consists of multiple skip convolutional gated recurrent unit blocks. The hidden state of the corresponding layer of the encoder is introduced into the decoder through skip connections.

8. The intelligent nearshore wave forecasting method suitable for digital simulation according to claim 7, characterized in that, The U-ConvGRU branch introduces a multi-scale wavelet transform layer to decompose the input boundary water level sequence into high-frequency detail signals and low-frequency trend signals, which are then sent to different channels. A gated routing weight mechanism is used at the neuron connection to dynamically trigger the hidden layer backtracking based on the real-time prediction error rate of the streaming data, and to perform local weight correction on time steps where the accumulated error exceeds the error threshold.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores program instructions, which, when executed in a computer, are used to perform the nearshore wave intelligent forecasting method suitable for digital simulation as described in any one of claims 1-8.

10. A nearshore wave intelligent forecasting system suitable for digital simulation, characterized in that, The system comprises the computer-readable storage medium of claim 9, wherein the system is a computer, the computer-readable storage medium is disposed within the system, and the system is provided with a microprocessor that executes program instructions stored in the computer-readable storage medium.

Images