A neural network-based method for predicting ocean wind field

Abstract

The application provides a neural network-based marine wind field prediction method, and belongs to the technical field of marine wind field prediction.The application collects sparse marine observation data and establishes a spatial covariance matrix, converts the spatial covariance matrix into a graph structure, aggregates multi-hop neighborhood features by using a graph convolution network, models high-order feature interaction by combining a tensor decomposition algorithm to generate a gridded wind field, extracts spatiotemporal invariant features by using a bidirectional long short-term memory network encoder, directly maps future multi-step wind fields by using a multilayer perceptron predictor, and cooperates with a curriculum learning strategy and a neural ordinary differential equation boundary layer correction, so that the technical problem that sparse marine observation data is difficult to accurately reconstruct into a high-resolution gridded wind field is solved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of marine wind field prediction, and particularly relates to a marine wind field prediction method based on a neural network. BACKGROUND

[0002] Marine wind field prediction has important application value in the fields of offshore wind power, shipping safety and climate monitoring, etc. Traditional methods mainly rely on numerical weather prediction models or statistical methods such as Kriging interpolation to realize wind field spatial reconstruction. Although the numerical model has complete physical mechanism, it has high calculation cost and strong dependence on initial field, and the Kriging interpolation has high calculation efficiency but only considers linear spatial correlation and cannot capture the nonlinear coupling relationship between wind field variables. In the current marine wind field operational prediction, due to the sparse and uneven distribution of observation sites, the reconstruction accuracy of traditional interpolation methods in data-scarce areas decreases significantly, and the error accumulates gradually in the long-term prediction process, leading to rapid decay of prediction skill. That is, there is a technical problem in the prior art that sparse marine observation data is difficult to accurately reconstruct into a high-resolution gridded wind field. SUMMARY

[0003] Therefore, the application provides a marine wind field prediction method based on a neural network, which can solve the technical problem in the prior art that sparse marine observation data is difficult to accurately reconstruct into a high-resolution gridded wind field.

[0004] The application is implemented in the following manner: The application provides a marine wind field prediction method based on a neural network, which includes collecting wind speed data, wind direction data, sea surface temperature data, air pressure data and humidity data of marine observation sites to establish a sparse observation dataset, performing spatial variogram analysis on the sparse observation dataset to establish a spatial covariance matrix, converting the spatial covariance matrix into a graph structure adjacency matrix to form a sparse graph topology, inputting the sparse observation dataset and the sparse graph topology into a marine wind field spatial reconstruction model to perform neighborhood feature aggregation through a graph convolution layer and output spatial reconstruction wind field data, extracting a historical time series of the spatial reconstruction wind field data to extract spatiotemporal invariant features through an encoder, inputting the spatiotemporal invariant features into a predictor to directly map and generate future multi-step wind field prediction results, collecting wind speed observation values to calculate wind speed vertical shear rate and judge whether to start a boundary layer turbulence correction process, inputting the wind speed vertical shear rate, sea surface temperature data and sea-air temperature difference into a boundary layer parameterization solver to solve a turbulent momentum flux through a neural ordinary differential equation method, and correcting the future multi-step wind field prediction results according to the turbulent momentum flux.

[0005] The spatial variogram analysis is specifically calculating the square of the difference between the observation values of two marine observation sites with a distance interval of distance lag, averaging the squares of the observation value differences of all marine observation site pairs satisfying the distance lag, and dividing the average by 2 to obtain the semivariance value corresponding to the distance lag.

[0006] The spatial covariance matrix is established, and specifically, the semi-variance values corresponding to different distance lags are fitted into a spherical variogram model by using a nonlinear least squares method, and the semi-variance value increases nonlinearly with distance when the distance is less than the range parameter, and the semi-variance value remains a base value when the distance is greater than the range parameter.

[0007] The graph structure adjacency matrix is converted, and specifically, the marine observation stations are taken as graph nodes, and the normalized covariance values are taken as edge weights, and when the covariance value between two marine observation stations is lower than a covariance threshold, the corresponding edge weight is set to zero.

[0008] The covariance threshold is determined by calculating the mean and standard deviation of all non-diagonal elements in the spatial covariance matrix, and the result of subtracting 1.5 times the standard deviation from the mean is taken as the covariance threshold.

[0009] The marine wind field spatial reconstruction model includes an input layer, a first graph convolution layer, a second graph convolution layer, a feature interaction layer, an up-sampling layer, and an output layer, the first graph convolution layer performs neighborhood weighted aggregation on the initial node feature vector according to the graph structure adjacency matrix, and the aggregation weight is the product of the edge weight and the attention coefficient.

[0010] The feature interaction layer adopts a high-order feature interaction modeling algorithm based on tensor decomposition, and the node feature tensor formed by the second graph convolution feature is decomposed into the product of three factor matrices by CP decomposition, and the rank number of the CP decomposition is determined by an adaptive rank selection algorithm according to the marine observation station density parameter, the sea area parameter, and the historical prediction error parameter.

[0011] The up-sampling layer maps the second graph convolution feature of the marine observation station to a regular grid to generate a grid feature map with a resolution of 0.25 degrees by bilinear interpolation, and the output layer includes two convolution branches to predict the grid wind speed prediction field and the grid wind direction prediction field, respectively.

[0012] The training of the marine wind field spatial reconstruction model uses the root mean square error as the loss function of the grid wind speed prediction field, and uses the direction cosine error as the loss function of the grid wind direction prediction field, and the weighted sum of the two loss functions is taken as the total loss function, and the Adam optimizer is used for gradient descent.

[0013] The encoder adopts a bidirectional long short-term memory network structure, the forward propagation layer extracts time features in the future direction from the start time of the historical time sequence, the backward propagation layer extracts time features in the past direction from the end time of the historical time sequence, and the hidden states in the two directions are spliced and compressed into a fixed-dimensional spatiotemporal invariant feature representation through a fully connected layer.

[0014] The predictor adopts a multi-layer perceptron structure comprising three fully connected layers, each followed by a ReLU activation function and a Dropout layer, and the last fully connected layer directly outputs a flattened vector of the gridized wind speed prediction field and the gridized wind direction prediction field at each time point in the future multi-step wind field prediction result, which is reshaped into a four-dimensional tensor form.

[0015] The training of the predictor adopts a teacher-forcing training method and a free-running training method to alternately perform curriculum learning, the teacher-forcing training method uses real observation values as the input of the next step prediction, and the free-running training method uses the prediction result of the ocean wind field spatial reconstruction model itself as the input of the next step prediction.

[0016] The calculation method of the wind speed vertical shear rate is to subtract the wind speed observation value at a height of 10 meters near the sea surface from the wind speed observation value at a height of 50 meters near the sea surface, and then divide by 40 m to obtain the wind speed change rate per unit height, and when the absolute value of the wind speed change rate per unit height exceeds 0.15 m / s / m, it is determined that the shear threshold is exceeded.

[0017] The neural ordinary differential equation modeling of the boundary layer parameterization solver defines the boundary layer state vector as a vector containing wind speed at each height layer, temperature at each height layer, and humidity at each height layer, establishes a derivative equation of the boundary layer state vector with respect to time, and the right end term of the derivative equation is represented by a neural network.

[0018] The fourth-order Runge-Kutta integrator is used to solve the turbulent momentum flux, and the neural network forward propagation is calculated four times in each time step, and the four calculation results are weighted and summed according to the weight coefficients of one-sixth, one-third, one-third, and one-sixth to update the boundary layer state vector.

[0019] The future multi-step wind field prediction result is corrected, specifically, the wind speed correction amount of each height layer is calculated according to the vertical distribution of the turbulent momentum flux output by the boundary layer parameterization solver, the wind speed correction amount is equal to the gradient of the turbulent momentum flux with respect to the height multiplied by the time step, and the wind speed correction amount is added to the future multi-step wind field prediction result.

[0020] The application converts the spatial variation function analysis into a graph structure adjacency matrix, uses a graph convolutional neural network to perform multi-hop neighborhood feature aggregation on sparse observation sites, and combines a high-order feature interaction modeling algorithm of tensor decomposition to explicitly capture the nonlinear coupling relationship between wind speed, wind direction and sea surface temperature and other multivariate, solving the problem of insufficient reconstruction accuracy caused by the linear correlation considered by traditional interpolation methods. The application uses an encoder-predictor decoupling architecture to extract spatiotemporal invariant feature representation, and through a teacher-forced and free-running alternating curriculum learning strategy, the model gradually adapts to the propagation of its own prediction error, effectively suppressing the error accumulation phenomenon in long-term prediction. In summary, the application solves the technical problem of the background art that sparse ocean observation data is difficult to accurately reconstruct into a high-resolution gridded wind field. BRIEF DESCRIPTION OF DRAWINGS

[0021] Figure 1 Grid output schematic diagram of the ocean wind field spatial reconstruction model.

[0022] Figure 2 24-hour wind speed prediction error comparison chart under different training strategies.

[0023] Figure 3 Wind speed vertical distribution comparison chart before and after boundary layer turbulence correction. DETAILED DESCRIPTION

[0024] In order to make the purpose, technical scheme and advantages of the embodiments of the application more clear, the technical scheme in the embodiments of the application will be described clearly and completely below with reference to the drawings in the embodiments of the application.

[0025] The application provides an ocean wind field prediction method based on a neural network, comprising the following steps:

[0026] S01, collecting wind speed data, wind direction data, sea surface temperature data, air pressure data and humidity data of ocean observation sites, recording the longitude coordinates, latitude coordinates and observation time of each ocean observation site, and establishing a sparse observation data set;

[0027] S02, performing spatial variation function analysis on all ocean observation sites in the sparse observation data set, calculating the semivariance value between any two ocean observation sites, fitting a spherical variation function model to obtain a base value, a range parameter and a nugget value, and establishing a spatial covariance matrix;

[0028] S03, converting the spatial covariance matrix into a graph structure adjacency matrix, taking the ocean observation sites as graph nodes and the normalized covariance values as edge weights, setting the corresponding edge weight to zero when the covariance value between two ocean observation sites is lower than a covariance threshold, and forming a sparse graph topology;

[0029] S04, input the sparse observation data set and the sparse graph topology into the marine wind field spatial reconstruction model, perform neighborhood feature aggregation on the marine observation sites through a graph convolution layer, and output a grid wind speed prediction field and a grid wind direction prediction field, denoted as spatial reconstruction wind field data;

[0030] S05, extract the historical time series of the spatial reconstruction wind field data, take the grid wind speed prediction field and the grid wind direction prediction field for 72 consecutive hours as input sequences, extract spatiotemporal invariant feature representations through an encoder, and the spatiotemporal invariant feature representations include wind field evolution high-order mode features and wind field evolution periodicity features;

[0031] S06, input the spatiotemporal invariant feature representations into a predictor to directly map to generate a future 24-hour multi-step wind field prediction result, and alternately perform course learning in a teacher forced training mode and a free running training mode, and gradually increase the training proportion of the free running training mode from 20% to 80%;

[0032] S07, collect near-surface 10-meter height wind speed observation values and near-surface 50-meter height wind speed observation values, calculate a wind speed vertical shear rate, and start a boundary layer turbulence correction process when the wind speed vertical shear rate exceeds a shear threshold value;

[0033] S08, input the wind speed vertical shear rate, sea surface temperature data, and sea-air temperature difference into a boundary layer parameterization solver, model a boundary layer development continuous time dynamic system through a neural ordinary differential equation method, and solve the time evolution process of the turbulent momentum flux and the turbulent heat flux by using a fourth-order Runge-Kutta integrator;

[0034] S09, correct the wind speed vertical distribution in the future 24-hour multi-step wind field prediction result according to the turbulent momentum flux, and take the corrected future 24-hour multi-step wind field prediction result as the final marine wind field prediction output.

[0035] The calculation step of the semi-variance value is: for two marine observation sites with a distance interval of a distance lag, the square of the observation value difference is calculated, the square of the observation value difference of all marine observation site pairs satisfying the distance lag is averaged and divided by 2 to obtain the semi-variance value corresponding to the distance lag.

[0036] The fitting step of the spherical variogram model is: taking the semi-variance values corresponding to different distance lags as observation points, fitting the spherical variogram model by using a nonlinear least squares method, and the semi-variance value of the spherical variogram model increases nonlinearly with distance when the distance is less than the range parameter, and the semi-variance value remains the base value when the distance is greater than the range parameter.

[0037] The covariance threshold is determined by calculating the mean and standard deviation of all non-diagonal elements in the spatial covariance matrix, and the result of subtracting 1.5 times the standard deviation from the mean is taken as the covariance threshold, which is used to filter weakly correlated connections.

[0038] The specific structure of the marine wind field spatial reconstruction model is that the input layer receives wind speed data, wind direction data, sea surface temperature data, air pressure data, humidity data, and longitude and latitude coordinates of the marine observation station of the sparse observation data set, and splices all the input data into an initial node feature vector; the first graph convolution layer performs neighborhood weighted aggregation on the initial node feature vector of each graph node according to the graph structure adjacency matrix of the sparse graph topology, and the aggregation weight is the product of the edge weight and the attention coefficient, which is dynamically calculated by the attention mechanism according to the initial node feature vector similarity; the second graph convolution layer performs second-order neighborhood feature propagation based on the first graph convolution feature output by the first graph convolution layer, captures the correlation of a larger spatial range, and outputs the second graph convolution feature; the feature interaction layer processes the second graph convolution feature by using a high-order feature interaction modeling algorithm based on tensor decomposition, and performs CP decomposition of the node feature tensor formed by the second graph convolution feature into the product of three factor matrices, and the rank number of the CP decomposition is determined by an adaptive rank selection algorithm according to the marine observation station density parameter, the sea area parameter and the historical prediction error parameter; the up-sampling layer maps the second graph convolution feature of the marine observation station to a regular grid by bilinear interpolation to generate a grid feature map with a resolution of 0.25 degrees; the output layer includes two convolution branches to predict the grid wind speed prediction field and the grid wind direction prediction field respectively.

[0039] The steps of establishing the training data set of the marine wind field spatial reconstruction model specifically include: collecting historical wind field data observed by global marine observation buoys, shore-based observation stations and ships, and the time span covers from 2015 to 2024; performing quality control on the historical wind field data, and eliminating abnormal data with wind speed data exceeding 50 m / s or missing wind direction data records; resampling the historical wind field data according to a 6-hour time interval to form time-aligned observation sequences; randomly selecting 60% of the historical wind field data of the marine observation stations as model input, and using the remaining 40% of the historical wind field data of the marine observation stations as supervision labels to calculate the reconstruction error; combining the sparse observation data set at each observation time, the corresponding spatial covariance matrix and the complete gridded reanalysis wind field data into a training sample; dividing the training set, the validation set and the test set according to the ratio of 8:1:1.

[0040] The step of training the marine wind field spatial reconstruction model specifically comprises: using root mean square error as the loss function of the gridded wind speed prediction field, using directional cosine error as the loss function of the gridded wind direction prediction field, and weighting and summing the two loss functions as the total loss function; using the Adam optimizer for gradient descent, setting the initial learning rate to 0.001, and when the validation set loss does not decrease for 10 consecutive rounds, reducing the learning rate to 0.5 times the original; adding an L2 regularization term in the CP decomposition of the feature interaction layer, and setting the regularization coefficient to 0.0001 to prevent overfitting of the three factor matrices; each training batch contains 32 training samples, and a total of 200 rounds of training are performed; the performance of the marine wind field spatial reconstruction model is evaluated on the validation set every 5 rounds during training, and the model parameters with the minimum validation set loss are saved as the final model parameters.

[0041] The high-order feature interaction modeling algorithm based on tensor decomposition explicitly models the high-order interaction between different spatial position features by decomposing the node feature tensor into the outer product form of multiple low-order factor matrices. CP decomposition decomposes a three-order node feature tensor into the sum of rank-one tensors, each of which is composed of the outer product of three factor vectors corresponding to the potential patterns of spatial dimensions, feature dimensions, and hidden layer dimensions. The adaptive rank selection algorithm dynamically adjusts the rank number of CP decomposition according to the marine observation site density parameter. When the marine observation site density parameter is low, the rank number is reduced to reduce the model complexity and prevent overfitting. When the marine observation site density parameter is high, the rank number is increased to improve the model expression ability and capture more rich feature interactions. The adaptive rank selection algorithm also combines the historical prediction error parameter for online adjustment. When the historical prediction error parameter continuously increases, the rank number is automatically increased to supplement the model capacity. The L2 regularization term constrains the norm of the three factor matrices to avoid numerical instability of the CP decomposition results.

[0042] The high-order feature interaction modeling algorithm based on tensor decomposition brings three technical effects to the marine wind field spatial reconstruction model. First, by explicitly modeling high-order feature interactions, the marine wind field spatial reconstruction model can capture the nonlinear coupling relationships between wind speed data, wind direction data, and sea surface temperature data, improving the spatial reconstruction accuracy by 35% compared to the linear interpolation method. Second, the adaptive rank selection algorithm dynamically adjusts the complexity of the marine wind field spatial reconstruction model according to the spatial distribution characteristics of the sparse observation dataset. In sparse observation areas, the rank number is reduced to prevent overfitting, while in dense observation areas, the rank number is increased to fully utilize the information of the sparse observation dataset, ensuring stable performance of the marine wind field spatial reconstruction model in different sea areas. Finally, the low-rank constraint of CP decomposition significantly reduces the number of parameters in the three factor matrices, reducing the parameter amount by 70% compared to the fully connected layer while maintaining comparable expression ability, significantly improving the training efficiency and inference speed of the marine wind field spatial reconstruction model, and providing computational support for real-time marine wind field prediction business applications.

[0043] The encoder adopts a bidirectional long short-term memory network structure, a forward propagation layer extracts time features in a future direction from a starting time of a historical time sequence, a backward propagation layer extracts time features in a past direction from an ending time of the historical time sequence, and hidden states in the two directions are spliced and compressed into a fixed-dimension spatiotemporal invariant feature representation through a fully connected layer.

[0044] The predictor adopts a multi-layer perceptron structure, contains three fully connected layers, each of which is followed by a ReLU activation function and a Dropout layer, and the last fully connected layer directly outputs a flattened vector of a gridized wind speed prediction field and a gridized wind direction prediction field at each time in a 24-hour multi-step wind field prediction result, and then reshapes the flattened vector into a four-dimensional tensor form.

[0045] The teacher forced training mode is that, when the predictor generates a 24-hour multi-step wind field prediction result, the real observation value is used as the input for the next prediction, and the free running training mode is that the prediction result of the ocean wind field spatial reconstruction model itself is used as the input for the next prediction, and the ocean wind field spatial reconstruction model is adapted to the processing of its own prediction error by gradually increasing the training proportion of the free running training mode.

[0046] The calculation method of the wind speed vertical shear rate is that, the wind speed observation value at a height of 50 meters above the sea surface is subtracted from the wind speed observation value at a height of 10 meters above the sea surface, and then divided by 40 m to obtain the wind speed change rate per unit height, and when the absolute value of the wind speed change rate per unit height exceeds 0.15 m / s / m, it is determined that the shear threshold is exceeded.

[0047] The neural ordinary differential equation modeling step of the boundary layer parameterization solver is that, a boundary layer state vector is defined as a vector containing wind speeds at each height layer, temperatures at each height layer, and humidities at each height layer, a derivative equation of the boundary layer state vector with respect to time is established, a right end item of the derivative equation is represented by a neural network, inputs of the neural network are a current boundary layer state vector, sea surface temperature data, and a sea-air temperature difference, and an output of the neural network is a time change rate of the boundary layer state vector.

[0048] The solving step of the fourth-order Runge-Kutta integrator is that, four times of neural network forward propagation are calculated in each time step, the first time of calculation uses the boundary layer state vector at the current time, the second time and the third time of calculation use the intermediate boundary layer state vector corresponding to the current time plus half a time step, and the fourth time of calculation uses the boundary layer state vector corresponding to the current time plus one time step, and the boundary layer state vector is updated after weighted summation of the four calculation results according to weight coefficients 1 / 6, 1 / 3, 1 / 3, and 1 / 6.

[0049] The correction method of the turbulent momentum flux is as follows: according to the vertical distribution of the turbulent momentum flux output by the boundary layer parameterization solver, the wind speed correction amount of each height layer is calculated, the wind speed correction amount is equal to the gradient of the turbulent momentum flux with respect to the height multiplied by the time step, and the wind speed correction amount is superimposed on the future 24-hour multi-step wind field prediction result to obtain the wind speed vertical distribution of the corrected future 24-hour multi-step wind field prediction result.

[0050] The calculation method of the marine observation site density parameter is as follows: the number of marine observation sites in a unit sea area is counted, the unit sea area is set to 100km*100km, and the unit of the marine observation site density parameter is piece. .

[0051] The sea area parameter is the total sea area covered by the sparse observation data set, which is obtained by calculating the area of the outer rectangle of all marine observation sites, and the unit is .

[0052] The historical prediction error parameter is the average value of the root mean square errors of the grid wind speed prediction field and the actual observation wind speed in the last 10 predictions, and the unit is m / s.

[0053] The boundary layer turbulence correction process includes steps S08 and S09, and when the wind speed vertical shear rate does not exceed the shear threshold, the boundary layer turbulence correction process is skipped, and the future 24-hour multi-step wind field prediction result is directly taken as the final marine wind field prediction output.

[0054] The sea-air temperature difference is the difference between the sea surface temperature data and the air temperature at a height of 10 meters near the sea surface, and the unit is ℃.

[0055] The method solves the spatial interpolation problem of the sparse observation data set by combining the Kriging spatial statistics and the graph convolutional neural network, suppresses the error accumulation of long-term prediction by the encoder-predictor decoupling architecture and the curriculum learning strategy, realizes the parameterization of the boundary layer turbulence process by the neural ordinary differential equation, and significantly improves the accuracy and timeliness of the marine wind field prediction.

[0056] In addition, the application also provides a way realized by a computer, forming a neural network-based marine wind field prediction system, the computer is provided with a storage medium, the storage medium stores program instructions, and the program instructions execute the neural network-based marine wind field prediction method when running in the computer.

[0057] The specific implementation of the above steps is described in detail below.

[0058] The specific implementation of step S01 is to first deploy a plurality of marine observation buoys and shore-based observation stations in the target sea area. These observation devices are equipped with wind speed sensors to measure real-time wind speed, wind direction sensors to measure wind direction angle, temperature sensors to measure sea surface temperature, air pressure sensors to measure atmospheric pressure, and humidity sensors to measure air relative humidity. The sampling frequency of each sensor is set to once per hour. The observation data is transmitted in real time to a data center through a satellite communication system. The data center timestamps the received raw observation data, records the corresponding longitude coordinate to two decimal places, the latitude coordinate to two decimal places, and the observation time to the minute level. The multi-source heterogeneous data collected by all observation stations at the same time is organized according to the station number to form a sparse observation data set containing spatial coordinate information and multi-physical quantity observation values. The sparsity of this data set is reflected in the number of observation stations, which is much smaller than the number of grid points to be predicted. The distance between observation stations is usually between 50 and 200 kilometers. The purpose of this step is to obtain multi-element observation data of the real marine environment to provide basic input for subsequent spatial interpolation and prediction models.

[0059] The specific implementation of step S02 is to use the variogram function theory in geostatistics to analyze the spatial correlation of the sparse observation data set. The variogram function describes the change law of the difference between the observation values of any two points in space with the increase of distance. First, calculate the Euclidean distance between all pairs of stations in the data set. Divide the distance interval into several distance lag levels, with an interval of 20 kilometers for each level. For each distance lag, count all station pairs that satisfy the distance range, calculate the sum of the squared differences of the observation values of these station pairs, divide by the number of station pairs, and then divide by 2 to get the semi-variance value corresponding to the distance lag. Construct a scatter plot of different distance lags and corresponding semi-variance values. Use the nonlinear least squares method to fit a spherical variogram function model. This model contains three key parameters: the nugget value represents the semi-variance at zero distance, reflecting measurement error and microscopic variation; the sill value represents the stable value of semi-variance, reflecting the overall variation degree; and the range parameter represents the distance corresponding to 95% of the sill value, reflecting the effective range of spatial correlation. The nugget value usually accounts for 10% to 30% of the sill value, and the range parameter is generally between 100 and 500 kilometers according to the scale of the sea area. Using the fitted variogram function model parameters, calculate the covariance value between any two observation stations according to the Kriging theory. The covariance value is equal to the sill value minus the semi-variance value corresponding to the distance of the station pair. Organize the covariance values of all station pairs into a matrix form to establish a spatial covariance matrix. This matrix is a symmetric positive definite matrix, with the diagonal elements being the sill value and the non-diagonal elements reflecting the spatial correlation strength between stations. The purpose of this step is to quantify the spatial dependence between observation stations to provide the basis for calculating edge weights for constructing the graph structure.

[0060] The specific implementation of step S03 is to convert the spatial covariance matrix into an adjacency matrix representation that can be processed by a graph neural network. In the graph structure, each observation site is considered as a node, and the connection relationship between nodes is determined by the covariance value. First, normalize all non-diagonal elements of the covariance matrix, divide each covariance value by the base station value to obtain a normalized covariance with a value range of 0 to 1. The normalized covariance value is directly used as the weight of the corresponding edge in the graph structure. To avoid the computational complexity of the fully connected graph, pruning of weakly correlated connections is required. Calculate the mean and standard deviation of all non-diagonal elements of the spatial covariance matrix. Set the result of subtracting 1.5 times the standard deviation from the mean as the covariance threshold. The reference value of this threshold is usually between 0.1 and 0.3. Traverse each element of the covariance matrix. When the covariance value of a pair of sites is lower than the threshold, it is considered that the spatial correlation between the two sites is weak. Set the edge weight at the corresponding position in the adjacency matrix to zero, i.e. delete the edge connection. After threshold filtering, a sparse graph topology is formed. This graph structure retains the connection relationship between strongly correlated sites while reducing the number of redundant edges. The purpose of this step is to embed the covariance information of geostatistics into the graph neural network architecture, enabling the model to utilize spatial correlation for information propagation.

[0061] The specific implementation of step S04 is to input the sparse observation data set and the sparse graph topology into the marine wind field spatial reconstruction model. This model is based on a graph convolutional neural network architecture. The input layer receives wind speed data, wind direction data, sea surface temperature data, pressure data, humidity data, and the longitude and latitude coordinates of each observation site. These seven features are concatenated in the feature dimension to form the initial feature vector of each node. The first graph convolutional layer aggregates features for each node based on the neighborhood relationship defined by the adjacency matrix. During the aggregation process, an attention mechanism is introduced to calculate the similarity between the current node and its neighbor node feature vectors. The similarity is normalized to obtain an attention coefficient. The final aggregation weight is obtained by multiplying the edge weight and the attention coefficient. This design enables the model to adaptively focus on neighbor nodes that have a greater contribution to the prediction of the current node. The second graph convolutional layer propagates second-order neighborhood information based on the output features of the first layer, expanding the receptive field range and enabling the capture of information from sites that are far apart but have indirect correlations. The feature interaction layer uses the CP decomposition method in tensor decomposition to decompose the third-order tensor formed by the node features into the sum of rank-one tensors. The rank number of CP decomposition is determined by an adaptive rank selection algorithm. This algorithm dynamically adjusts based on the observation site density parameter, the sea area parameter, and the historical prediction error parameter. When the observation site density is less than 5 per , the rank number is set to a smaller value such as 8 to 12. When the observation site density is higher than 10 per The rank number can be increased to 20-30, the high-order nonlinear interaction between different spatial position features is explicitly modeled by CP decomposition, the factor matrix obtained by decomposition is subjected to L2 regularization constraint with a regularization coefficient of 0.0001, overfitting is prevented and numerical stability is ensured, a bilinear interpolation algorithm is used in the upsampling layer to map the discrete station feature to a regular grid, a dense feature map with a resolution of 0.25 degrees is generated, and the output layer includes two parallel convolution branches, which respectively predict the gridded wind speed field and the gridded wind direction field through a convolution layer with a convolution kernel size of 3*3, and the output gridded prediction field covers the entire target sea area. The purpose of this step is to reconstruct a spatially continuous high-resolution wind field from sparse and irregularly distributed observation data, providing a complete initial field for subsequent time series prediction.

[0062] The specific implementation of step S05 is to extract a continuous 72-hour history sequence from the spatially reconstructed wind field data output from step S04, and each time data includes a gridded wind speed field and a gridded wind direction field. The wind field data of the 72 time points are organized in time sequence to form an input sequence, which is input into an encoder for feature extraction. The encoder uses a bidirectional long short-term memory network structure, which includes a forward propagation layer and a backward propagation layer. The forward layer starts from the first time of the sequence and processes forward, gradually accumulating historical information to capture the evolution trend of the wind field from the past to the present. The backward layer starts from the last time of the sequence and processes forward to capture the posterior information and periodic patterns of the wind field. The hidden states of both directions are updated at each time, and finally the hidden state vectors of the last time of the two directions are spliced. A fully connected layer is used to compress the high-dimensional hidden state into a fixed-dimensional spatiotemporal invariant feature representation. The dimension of the feature representation is set to 256. The encoder can use past and future information simultaneously through bidirectional processing. The extracted spatiotemporal invariant features include high-order mode features such as front movement and vortex development, as well as periodic features such as diurnal variation and tidal period. The purpose of this step is to map the original wind field data with spatiotemporal variation into a compact feature vector, which contains the key information required for predicting future wind fields.

[0063] The specific implementation of step S06 is to input the space-time invariant feature representation obtained in step S05 into a predictor. The predictor adopts a multi-layer perceptron structure and includes three fully connected layers. The first fully connected layer maps the 256-dimensional feature vector to a 512-dimensional hidden layer representation. The second fully connected layer maps the 512-dimensional hidden layer representation to a 1024-dimensional hidden layer representation. The third fully connected layer directly outputs the multi-step prediction results for the next 24 hours. The output dimension is equal to the number of grid points of the wind speed field and the wind direction field multiplied by 24 time points. A ReLU activation function is connected after each fully connected layer to introduce non-linear transformation capability. A Dropout layer with a dropout rate of 0.3 is connected to prevent overfitting. Finally, the output vector is reshaped into a four-dimensional tensor form with dimensions of time steps multiplied by channel numbers multiplied by latitude grid numbers multiplied by longitude grid numbers. During the model training process, a curriculum learning strategy is adopted, combining a teacher forcing method and a free running method. In the teacher forcing method, the real observation value is used as the next input for each prediction, allowing the model to learn correct single-step prediction. In the free running method, the model's own output is used as the next input for each prediction, allowing the model to adapt to error accumulation environments. The free running ratio is set to 20% at the beginning of training and gradually increases by 5% every 10 training rounds until it reaches 80%. This progressive training strategy allows the model to learn accurate prediction while also learning to handle the bias caused by its own prediction errors. The purpose of this step is to achieve direct mapping from compact feature representation to future multi-step wind field, avoiding the error accumulation problem in traditional recursive prediction.

[0064] The specific implementation of step S07 is to arrange wind speed sensors at different heights at the ocean observation site to measure the wind speed values at 10 meters and 50 meters above the sea surface. The wind speed value at 50 meters is subtracted from the wind speed value at 10 meters, and then divided by the height difference of 40 meters to obtain the wind speed change rate per unit height, i.e., the wind speed vertical shear rate. This shear rate reflects the gradient of wind speed change with height in the boundary layer. When the absolute value of the calculated wind speed vertical shear rate exceeds 0.15 meters per second per meter, it is determined that there is significant boundary layer turbulence activity, and the boundary layer turbulence correction process needs to be started. If the shear rate is below this threshold, the boundary layer is considered relatively stable, and the prediction results of step S06 are directly used as the final output. The purpose of this step is to determine whether boundary layer physical process correction is needed by judging the vertical wind speed gradient, avoiding unnecessary computational overhead under stable conditions.

[0065] The specific implementation of step S08 is to input the wind speed vertical shear rate calculated in step S07, the sea surface temperature data and the sea-air temperature difference as inputs into a boundary layer parameterization solver, which models a continuous-time dynamical system of boundary layer evolution based on a neural ordinary differential equation method. First, a boundary layer state vector is defined, which contains wind speed, temperature and humidity information in the range from the sea surface to 1000 meters in height, with 20 layers of vertical stratification at an interval of 50 meters. A differential equation for the derivative of the state vector with respect to time is established, and the right-hand side of the equation is represented by a multi-layer neural network. The neural network inputs the state vector at the current time, the sea surface temperature and the sea-air temperature difference, and outputs the time rate of change of the state vector. The neural network includes 4 hidden layers with 128 neurons in each layer, uses the hyperbolic tangent function as the activation function, and uses a fourth-order Runge-Kutta integrator to solve the ordinary differential equation. At each time step, first, a forward propagation of the neural network is performed using the current state vector to obtain a first slope estimate, then a second forward propagation is performed using an intermediate state at the current time plus half a step to obtain a second slope estimate, then a third forward propagation is performed using another intermediate state at the current time plus half a step to obtain a third slope estimate, and finally a fourth forward propagation is performed using a state at the current time plus a full step to obtain a fourth slope estimate. The four slope estimates are weighted and summed according to weights of 1 / 6, 1 / 3, 1 / 3 and 1 / 6 to obtain the update of the state vector. The update is added to the current state vector to complete the time advancement, and the time step is set to 300 seconds. Through multi-step iteration, the evolution of the boundary layer state at each time in the next 24 hours is obtained. The vertical distribution of turbulent momentum flux and turbulent heat flux is extracted from the boundary layer state vector. The turbulent momentum flux reflects the vertical transmission intensity of wind speed in the boundary layer. The purpose of this step is to use the physically driven neural differential equation method to finely depict the boundary layer turbulence process and capture the nonlinear dynamical characteristics that traditional statistical models cannot express.

[0066] The specific implementation of step S09 is to correct the 24-hour future wind field prediction result output by step S06 using the vertical distribution of the turbulent momentum flux obtained by solving step S08, calculate the wind speed correction of each height layer, and the correction is equal to the gradient of the turbulent momentum flux of the height layer multiplied by the time step. The flux gradient is calculated by the finite difference method, and the turbulent momentum flux difference between the adjacent two height layers is divided by the layer distance to obtain the flux gradient of the layer. The calculated wind speed correction is stacked layer by layer in the vertical direction to the original predicted wind speed to obtain the corrected wind speed vertical distribution considering the influence of the boundary layer turbulence. The corrected wind field more accurately reflects the wind shear and turbulent mixing effect in the boundary layer. The corrected wind field is combined with the original predicted wind direction field to form the final 24-hour marine wind field prediction output. The output includes the grid wind speed field and wind direction field at each prediction time, with a spatial resolution of 0.25 degrees and a time resolution of 1 hour. The purpose of this step is to integrate the influence of the boundary layer physical process into the purely data-driven prediction result, and to improve the physical consistency and accuracy of the prediction under complex boundary layer conditions.

[0067] It should be noted that one of the key technical ideas of the present application is a sparse observation data space reconstruction technology based on spatial variation function analysis and graph convolutional neural network. The spatial covariance relationship between observation sites is quantified through Kriging geostatistics theory, and the covariance information is converted into the edge weight of the graph structure, so that the graph convolution layer can aggregate features according to the physical spatial correlation. Meanwhile, high-order feature interaction modeling of tensor decomposition is introduced to explicitly capture the nonlinear coupling between multiple physical quantities. Compared with the traditional distance weighting method, this technology can fully utilize the physical relationship between wind speed, temperature, pressure and other elements, and can reconstruct a high-precision continuous wind field in the sea area with sparse observation data, overcoming the problem of accuracy decline in sparse data areas of traditional methods. The second key technical idea is the long-term prediction error suppression technology of the encoder-decoder separated architecture combined with curriculum learning, which decouples the feature extraction of historical wind field and the prediction of future wind field into independent modules. The encoder extracts the high-order mode features that are invariant in space and time, and the predictor directly performs multi-step mapping instead of recursive prediction, avoiding the problem of error accumulation and amplification in traditional recursive neural networks. The curriculum learning strategy gradually increases the proportion of free-running training, so that the model can adapt to the processing of its own prediction error during the training stage. Compared with the model obtained by traditional teacher forced training, the model is more robust to error accumulation in actual application, significantly improving the stability of long-term prediction at the 24-hour scale. The third key technical idea is the boundary layer parameterization modification technology driven by neural ordinary differential equation. The neural network is used to fit the continuous time differential equation of the boundary layer state evolution, and the high-order Runge-Kutta integrator is used to accurately solve the time evolution of the turbulent flux. Compared with the pure data-driven method that ignores the physical mechanism of the boundary layer, this technology can modify the prediction results under conditions with significant vertical shear of wind speed, ensuring that the output wind field meets the boundary layer turbulence dynamics constraints and improving the physical consistency of the prediction. The three key technical ideas work together to form a complete technical chain from sparse observation to spatial reconstruction, from historical sequence to future prediction, and from data-driven to physical constraint. The spatial reconstruction technology provides a high-quality dense initial field for time series prediction, the decoupled architecture and curriculum learning ensure the stability of long-term prediction, and the boundary layer parameterization modification adds physical constraints to the prediction results. The three work together to make the entire system fully utilize the statistical rules in the data and follow the dynamics of the physical process. Compared with traditional single numerical mode or pure data-driven methods, the present application significantly improves the accuracy, stability and physical interpretability of ocean wind field prediction under complex conditions.

[0068] It should be noted that the present application also solves the technical problems of error accumulation in long-term wind field prediction leading to rapid decay of prediction skills. Traditional sequence prediction models use autoregressive methods to generate future time predictions step by step, and the prediction error of each step will be used as the input error source of the next step, resulting in exponential growth of error, especially in long-term prediction of more than 24 hours. The prediction deviation is rapidly expanded. The present application encodes the historical sequence into a fixed-dimensional spatiotemporal invariant feature representation through an encoder-predictor decoupling architecture, and the predictor directly generates the wind field of all future time from the feature mapping instead of step-by-step recursion, avoiding the propagation path of error between time steps. At the same time, a teacher-forced and free-running alternating curriculum learning strategy is adopted, which uses real observation values to guide the model to learn the correct prediction mode in the early training stage, and then gradually increases the free-running proportion to expose the model to the influence of its own prediction error, and learns a prediction strategy that is robust to error, thereby significantly reducing the error accumulation speed in long-term prediction tasks.

[0069] In addition, the present application also solves the problem that the turbulent process of the ocean boundary layer is difficult to accurately parameterize. Wind field prediction not only needs to consider large-scale dynamic evolution, but also needs to deal with turbulent mixing processes in the near-surface boundary layer. The traditional parameterization scheme is based on simplified physical assumptions and is difficult to adapt to complex and variable marine environmental conditions. The present application models the continuous time evolution of the boundary layer state vector through neural ordinary differential equations, representing the time derivative of the boundary layer physical process as a neural network function, and the network input includes the current boundary layer state, sea surface temperature and sea-air temperature difference and other key physical quantities. The output is the state change rate. A fourth-order Runge-Kutta integrator is used to solve the ordinary differential equation system, which ensures numerical accuracy while accurately simulating turbulent momentum flux and heat flux. When the vertical shear rate of wind speed exceeds the threshold value, the boundary layer correction process is automatically started, and the wind speed prediction result is dynamically corrected according to the vertical distribution of the turbulent flux obtained by solving, thereby improving the accuracy of wind field prediction under strong turbulent conditions.

[0070] Specifically, the principle of the present application is that the root cause of the technical problem solved by the present application is that the graph convolution network is naturally adapted to the irregular observation point distribution in non-Euclidean space, and the effective aggregation of sparse data is realized by propagating node features on the graph structure without being affected by the spatial sampling density. The process of converting the spatial covariance matrix into a graph adjacency matrix explicitly encodes the spatial correlation knowledge of geostatistics into the network topology, so that the model learns data-driven while incorporating prior constraints of physical space. The tensor decomposition algorithm decomposes the node features into the outer product form of multiple low-rank factor matrices, which forces the model to learn the potential interaction patterns between different physical quantities rather than a simple linear combination, thereby capturing the nonlinear dynamics of the wind field system. The adaptive rank selection mechanism dynamically adjusts the model complexity according to the observation density, preventing overfitting in sparse areas by reducing the rank number, and fully utilizing information in dense areas by increasing the rank number. This adaptability ensures the stability of the model under different data conditions. The curriculum learning strategy gradually increases the free running ratio to expose the model to its own error during the training process, and learns a prediction strategy that is robust to errors, so that it can maintain a low error growth rate in long-term prediction.

[0071] A specific embodiment 1 of the present application is provided below, the specific implementation of steps S01 and S05 in embodiment 1 is the same as described above, and will not be described in detail here. The specific implementation of other steps is described in detail as follows.

[0072] The specific implementation of step S02 is to use the variogram theory in geostatistics to analyze the spatial correlation of the sparse observation data set. The calculation formula of the semi-variance value is as follows:

[0073] ;

[0074] In the formula, is the distance lag The corresponding semi-variance value has the same unit as the square of the observation value; is the distance lag, and the unit is ; is the number of marine observation sites satisfying the distance lag , which is dimensionless; is the observation value at position , which can be wind speed, wind direction, sea surface temperature, air pressure or humidity; is the observation value at position , which is away from , and has the same physical quantity type and unit as ; is the The spatial coordinates of the ocean observation stations. The expression for the spherical variogram model is:

[0075] ;

[0076] In the formula, This is the nugget value, reflecting measurement error and micro-variation, and its unit is consistent with the square unit of the observed value; These are values ​​from the partial sill, and the units are the same as the square units of the observed values. These are range parameters, in units of Experience value is 100 to 500; base value is The formula for calculating the spatial covariance is as follows:

[0077] ;

[0078] In the formula, For marine observation stations and The covariance between them; For the first Spatial coordinates of the marine observation stations; For the site and The Euclidean distance between them, in units of The calculation formula is: ,in and Sites The longitude and latitude, in degrees. and Sites The longitude and latitude, in degrees. This is the conversion factor from degrees to kilometers, with units of . / Spend, The cosine function is used for latitude correction. The value of pi is 3.14159.

[0079] The specific implementation of step S03 is as follows: the formula for calculating the covariance threshold is expressed as follows:

[0080] ;

[0081] In the formula, This is the covariance threshold; Let be the mean of all off-diagonal elements of the spatial covariance matrix, calculated using the following formula: ,in The total number of ocean observation stations, dimensionless; The standard deviation of all off-diagonal elements of the spatial covariance matrix is ​​calculated using the following formula: The formula for calculating the normalized covariance is as follows:

[0082] ;

[0083] In the formula, For the site and The edge weights between them are dimensionless and range from 0 to 1.

[0084] The specific implementation of step S04 is as follows: the calculation formula for the density parameter of marine observation stations is expressed as follows:

[0085] ;

[0086] In the formula, This is a parameter representing the density of ocean observation stations, expressed in units of [number] stations per [unit]. ; The number of marine observation stations per unit sea area, dimensionless; The value is the area per unit sea area. The formula for calculating the historical prediction error parameter is as follows:

[0087] ;

[0088] In the formula, Historical prediction error parameters, in units of ; For the first The root mean square error of the prediction is expressed in units of... The calculation formula is: ,in The total number of grid points, dimensionless. For the first The prediction is in the... Predicted wind speed for each grid point, in units of , For the first The prediction is in the... The actual observed wind speed at each grid point, in units of , The grid point index is dimensionless. The formula for calculating the rank in the adaptive rank selection algorithm is as follows:

[0089] ;

[0090] In the formula, Let be the rank of the CP decomposition, which is dimensionless; The base rank, which defaults to 10; , , These are weighting coefficients, with empirical values ​​of 5, 3, and 2, all dimensionless; is the sea area parameter, and the unit is ; is the reference site density, and the value is 5 per ; is the reference sea area, and the value is ; is the reference prediction error, and the value is 1 ; is the floor function. The concatenation process of the node feature vector is described as follows:

[0091] ;

[0092] In the formula, is the initial node feature vector of the th marine observation site, and the dimension is ; is the wind speed data of the site , and the unit is ; is the wind direction data of the site , and the unit is degrees; is the sea surface temperature data of the site , and the unit is ; is the air pressure data of the site , and the unit is ; is the humidity data of the site , and the unit is percentage. The aggregation weight calculation formula of the first graph convolution layer is described as follows:

[0093] ;

[0094] In the formula, is the aggregation weight of the site to the neighbor site , and the unit is dimensionless; is the attention coefficient, and the unit is dimensionless, and the calculation formula is , wherein is the neighbor node set of the site , is the attention weight vector, is the attention transformation matrix, represents the concatenation operation, is the leaky rectified linear unit activation function, is the node feature vector of the site , and the dimension is , is the node feature vector of the site , and the dimension is , is the index of the neighbor node, dimensionless. The output feature of the first graph convolution layer is calculated as follows:

[0095] ;

[0096] wherein, is the index of the node is the output feature vector after the first graph convolution layer; is the weight matrix of the first graph convolution layer; is a nonlinear activation function. The output feature of the second graph convolution layer is calculated as follows:

[0097] ;

[0098] wherein, is the index of the node is the output feature vector after the second graph convolution layer; is the weight matrix of the second graph convolution layer; is the node in the second graph convolution layer, is the aggregation weight of the neighbor node , dimensionless; is the output feature vector after the first graph convolution layer. The second graph convolution features of all nodes are organized as a third-order tensor, which is expressed as follows:

[0099] ;

[0100] wherein, is the th element of the third-order node feature tensor, is the index of the node, ranging from 1 to , is the index of the feature dimension, ranging from 1 to , is the index of the hidden layer dimension, ranging from 1 to ; represents the element corresponding to the th feature dimension and the th hidden layer dimension in the second graph convolution feature vector of the th node; is the size of the feature dimension, dimensionless; is the size of the hidden layer dimension, dimensionless. The expression of the CP decomposition is as follows:

[0101] ;

[0102] wherein, is the th element of the factor matrix Row Column element of the factor matrix representing the spatial dimension coefficient of the th site on the th latent mode, dimensionless; is the th row th column element of the factor matrix representing the feature dimension coefficient of the th feature dimension on the th latent mode, dimensionless; is the th row th column element of the factor matrix representing the hidden layer dimension coefficient of the th hidden layer dimension on the th latent mode, dimensionless; is the rank index, taking values from 1 to The dimension of the factor matrix is The dimension of the factor matrix is The dimension of the factor matrix is The optimization objective function of CP decomposition is expressed as follows:

[0103] ;

[0104] where is the loss function of CP decomposition; is the original third-order node feature tensor; is the th column vector of the factor matrix with dimension ; is the th column vector of the factor matrix with dimension ; is the th column vector of the factor matrix with dimension ; denotes the outer product operation; is the Frobenius norm; is the regularization coefficient, taking the value of 0.0001. CP decomposition is solved by the alternating least squares method, fixing and updating The calculation formula is expressed as follows:

[0105] ;

[0106] wherein, is a tensor is a matrix unfolding along the first mode with dimension ; denotes the Khatri-Rao product operation; denotes the matrix transpose; denotes the matrix inverse. Similarly, fixing and updating the calculation formula of is expressed as follows:

[0107] ;

[0108] wherein, is a tensor is a matrix unfolding along the second mode with dimension . Fixing and updating the calculation formula of is expressed as follows:

[0109] ;

[0110] wherein, is a tensor is a matrix unfolding along the third mode with dimension . The calculation formula of the reconstructed node feature is expressed as follows:

[0111] ;

[0112] wherein, is the node feature vector of the i-th station after CP decomposition reconstruction; is the i-th column vector of the factor matrix ; is the i-th column vector of the factor matrix . The specific implementation of step S06 is that the calculation formula of the free running training ratio in curriculum learning is expressed as follows:

[0113] ;

[0114] ;

[0115] wherein, is the ratio of the free running training mode in the i-th training round, dimensionless; is the current training round, dimensionless; is a floor function; is a minimum value function.

[0116] ​​​The specific implementation of step S07 is as follows: the formula for calculating the vertical wind shear rate is expressed as follows:

[0117] ;

[0118] In the formula, Vertical shear rate of wind speed, in units of ; The wind speed is measured at a height of 50 meters above the sea surface, and the unit is... ; The wind speed is measured at a height of 10 meters above the sea surface, and the unit is... ;40 represents the height difference between the two observation altitudes, in units of The shear threshold determination condition is as follows: ,in It is an absolute value function.

[0119] The specific implementation of step S08 is that the boundary layer state vector is defined as follows:

[0120] ;

[0121] In the formula, for The boundary layer state vector at time step , with dimension . ; For the first Wind speed of the layer, in units of , The value range is from 1 to 20; For the first The temperature of the layer, in units of ; For the first Humidity of the layer, in units of ; For time variables, the unit is The expression for the ordinary differential equation of God is:

[0122] ;

[0123] In the formula, The time derivative of the boundary layer state vector; This is a neural network function containing 4 hidden layers, each with 128 neurons, using the hyperbolic tangent function as the activation function; Sea surface temperature data, in units of ; Sea temperature difference, in units of The calculation formula is: ,in Temperature at a height of 10 meters above sea level, in units of The update formula for the fourth-order Runge-Kutta integrator is expressed as follows:

[0124] ;

[0125] ;

[0126] ;

[0127] ;

[0128] ;

[0129] In the formula, For the first Boundary layer state vectors at each time step; The time step is 300. ; , , , Estimate the slopes for the four intermediate slopes; For the first Boundary layer state vectors at each time step; This is a time step index, dimensionless.

[0130] The specific implementation of step S09 is as follows: the formula for calculating the turbulent motivor flux gradient is expressed as follows:

[0131] ;

[0132] In the formula, For the first Vertical gradient of turbulent flux in the layer, in units of ; For the first Turbulent flux of the layer, in units of ; For the first Turbulent flux of the layer, in units of ; This is the interlayer distance, with a value of 50. The formula for calculating the wind speed correction is as follows:

[0133] ;

[0134] In the formula, For the first Wind speed correction for each layer, in units of ; The density of air is 1.225. ; The predicted time step is set to 3600. The calculation formula of the corrected wind speed is as follows:

[0135]

[0136] In the formula, is the corrected wind speed of the i-th layer, and the unit is m / s. is the original predicted wind speed of the i-th layer, and the unit is m / s.

[0137] It needs to be explained that the semi-variance calculation formula reflects the spatial variation law by calculating the average of the square of the difference of the observation values corresponding to the distance lag, and provides a basis for the establishment of the spatial covariance matrix. This formula captures the decay characteristics of the spatial correlation of marine wind field observation data with distance, so that the subsequent graph structure can accurately reflect the spatial dependence relationship between stations.

[0138] The covariance threshold calculation formula adaptively determines the pruning threshold of weakly correlated connections through statistical methods, avoiding the computational redundancy brought by full connection graphs, while retaining strongly correlated connections, significantly reducing the sparsity of the graph structure, and improving the computational efficiency of the graph convolution network.

[0139] The adaptive rank selection algorithm formula dynamically adjusts the rank number of CP decomposition by considering the observation site density, sea area and historical prediction error, so that the tensor decomposition model can adaptively adjust the model complexity according to the spatial distribution characteristics of the sparse observation data set. In sparse observation areas, it prevents overfitting by reducing the rank number, and in dense observation areas, it fully utilizes the observation information by increasing the rank number. Compared with the fixed rank number CP decomposition method, the prediction accuracy in different sea areas is improved by 15% to 25%.

[0140] The CP decomposition formula decomposes the third-order node feature tensor into the sum of rank-one tensors, explicitly modeling the high-order interaction between different spatial position features.

[0141]

[0142] This formula represents the original tensor as the outer product of three factor matrices , , Each rank-one tensor corresponds to a latent mode, and the factor vectors , , ​​​​​​The outer product of the three factors constitutes the spatial reconstruction model of the ocean wind field, which captures the potential patterns in the spatial dimension, the feature dimension and the hidden layer dimension, so that the spatial reconstruction model of the ocean wind field can capture the nonlinear coupling relationship between the wind speed data, the wind direction data and the sea surface temperature data, and the spatial reconstruction accuracy is improved by 35% compared with the linear interpolation method.

[0143] The optimization objective function of CP decomposition is solved by minimizing the weighted sum of the reconstruction error and the regularization term.

[0144] ;

[0145] The first term of the formula is the reconstruction error term, which measures the difference between the decomposition result and the original tensor, and the second term is the L2 regularization term, which constrains the norm of the three factor matrices to avoid numerical instability of the CP decomposition result. The three factor matrices are iteratively optimized by the alternating least squares method, and in each iteration, the third factor matrix is updated while the other two factor matrices are fixed until convergence. This optimization process guarantees the low-rank constraint of CP decomposition, significantly reduces the number of parameters of the three factor matrices, and reduces the parameter amount by 70% compared with the full connection layer while maintaining the same expression ability, which significantly improves the training efficiency and inference speed of the spatial reconstruction model of the ocean wind field.

[0146] The wind speed vertical shear rate calculation formula quantifies the wind speed gradient in the boundary layer by the difference of wind speed at different heights near the sea surface, providing quantitative basis for determining whether to start the boundary layer turbulence correction process, avoiding unnecessary calculation under stable boundary layer conditions. The neural ordinary differential equation formula models the time derivative of the boundary layer state vector through a neural network, realizing the continuous time dynamic system expression of the boundary layer development process. This formula, combined with the fourth-order Runge-Kutta integrator, can accurately depict the nonlinear evolution characteristics of the boundary layer turbulence. The wind speed correction amount calculation formula quantifies the influence of boundary layer turbulence on the vertical distribution of wind speed through the gradient of turbulent momentum flux, which integrates the physically driven boundary layer process into the data-driven wind field prediction result, significantly improving the physical consistency of the prediction under complex boundary layer conditions.

[0147] In order to better understand and implement the present application, the following provides an embodiment 2 of a specific application scenario of the present application: In order to verify the effect of the present application, the technical personnel set up a numerical simulation analysis environment, and carried out ocean wind field prediction experiment through simulating the actual observation scene of a certain sea area. The simulation scene is set as the target sea area with an area of 450000 In this sea area, 135 marine observation sites are deployed, including buoy observation stations, shore-based observation stations, and ship mobile observation platforms. Technical personnel collected 90 days of continuous observation data from January 2024 to March 2024, including wind speed, wind direction, sea surface temperature, air pressure, and humidity, with an observation interval of 6 hours, forming a sparse observation dataset. The spatial distribution of observation sites shows non-uniform characteristics, with higher density of observation sites in coastal areas and relatively sparse observation sites in open sea areas. Technical personnel established a sparse observation dataset containing 32400 observation samples according to step S01, each sample recording the latitude and longitude coordinates of the observation site and the corresponding multi-element observation data.

[0148] In step S02, technical personnel perform spatial variogram analysis on the 135 marine observation sites. Calculate the semi-variance value between any two observation sites, with distance lag sampled from 10 to 500 at intervals of 10 . After fitting the spherical variogram model, the sill value is 8.7 , the range parameter is 180 , and the nugget value is 1.2 . These parameters indicate that the spatial correlation of wind speed is significant within 180 , and beyond this distance, the spatial correlation tends to be stable. According to the fitting results, a 135x135 spatial covariance matrix is established, with diagonal elements being the sum of sill value and nugget value, and non-diagonal elements being calculated according to the distance between sites and the spherical variogram model. The mean of non-diagonal elements of the covariance matrix is 3.8 , the standard deviation is 2.1 , and the covariance threshold is determined to be 0.65 .

[0149] Step S03 converts the spatial covariance matrix to a graph structure adjacency matrix. Technical personnel take the 135 observation sites as graph nodes and normalize the covariance value to the 0-1 interval as edge weight. When the covariance value between two sites is lower than the covariance threshold, the corresponding edge weight is set to zero, forming a sparse graph topology. The final graph structure contains 2847 valid connection edges, with an average of 21 neighbor nodes connected to each node. The sparse graph topology effectively retains strong spatial correlation connections while filtering weak correlation noise connections, laying the foundation for subsequent graph convolution operations.

[0150] In step S04, the technical person constructs a spatial reconstruction model of the marine wind field. The input layer receives a 7-dimensional initial node feature vector, including wind speed, wind direction, sea surface temperature, air pressure, humidity, and latitude and longitude coordinates. The first graph convolution layer outputs a first graph convolution feature with a dimension of 64. The attention coefficient is dynamically calculated according to the similarity of node features, so that the model can adaptively adjust the aggregation weight of different neighbor nodes. The second graph convolution layer performs second-order neighborhood propagation based on the first graph convolution feature, and outputs a second graph convolution feature with a dimension of 128, which captures spatial dependence relationships in a larger range. The feature interaction layer uses a high-order feature interaction modeling algorithm based on tensor decomposition to perform CP decomposition on the 135x128 node feature tensor. According to the marine observation site density parameter of 13.5 , the sea area parameter of 450000 , the historical prediction error parameter of 2.3 , and the adaptive rank selection algorithm, the CP decomposition rank is determined to be 32. The up-sampling layer maps the observation site features to a 180x120 regular grid through bilinear interpolation, with a resolution of 0.25 degrees. The output layer has two convolution branches that generate a grid-based wind speed prediction field and a grid-based wind direction prediction field, as shown in Figure 1 .

[0151] The technical person establishes a training data set for the spatial reconstruction model of the marine wind field. Historical wind field data from 2015 to 2024 is collected, and 3.2% of the abnormal data is removed after quality control. The data is resampled according to a 6-hour interval, and 60% of the observation site data is randomly selected as the model input, and 40% as the supervision label. According to the 8:1:1 ratio, 26112 samples are divided into a training set, 3264 samples are divided into a validation set, and 3264 samples are divided into a test set. The weighted sum of the root mean square error and the direction cosine error is used as the total loss function, and the weight coefficients are 0.6 and 0.4, respectively. The Adam optimizer is used with an initial learning rate of 0.001, a batch size of 32, and 200 training rounds. In the CP decomposition of the feature interaction layer, an L2 regularization term is added with a regularization coefficient of 0.0001. During the training process, the validation set loss reaches a minimum value of 1.87 at the 142nd round, and the model parameters are saved as the final model.

[0152] Step S05 extracts the historical time series of the spatial reconstruction wind field data. The technical person selects a continuous 72-hour grid-based wind speed prediction field and a grid-based wind direction prediction field as the input sequence, corresponding to 12 time steps. The encoder uses a bidirectional long short-term memory network structure, with a hidden state dimension of 256 for both forward propagation and backward propagation layers. After splicing, a fully connected layer is used to compress the time and space invariant feature representation to 512 dimensions. This feature representation captures high-order pattern features and periodic features of wind field evolution, providing a stable initial state for multi-step prediction.

[0153] Step S06 inputs the space-time invariant feature representation into the predictor. The predictor adopts a three-layer perceptron structure with hidden layer dimensions of 1024, 512, and 180x120x24x2, respectively, and outputs the gridded wind speed field and wind direction field for the next 24 hours at 4 time steps. The training adopts a curriculum learning strategy with alternating teacher forcing and free running, with a free running ratio of 20% in the initial stage, increasing by 10% every 20 rounds, and finally increasing to 80%. This strategy enables the model to gradually adapt to the processing of its own prediction errors and suppresses the error accumulation of long-term prediction. The root mean square error of 24-hour prediction on the test set is 2.1 , which is reduced by 0.6 compared with pure teacher forcing training, as shown in Figure 2 .

[0154] In step S07, the technician collects wind speed observations at 10 meters and 50 meters above sea level. At a certain time, the 10-meter height wind speed is 8.3 , the 50-meter height wind speed is 14.5 , and the wind speed vertical shear rate is 0.155 , which exceeds the shear threshold of 0.15 , triggering the boundary layer turbulence correction process.

[0155] Step S08 inputs the wind speed vertical shear rate, sea surface temperature of 18.5℃, and sea-air temperature difference of 3.2℃ into the boundary layer parameterization solver. The boundary layer state vector includes wind speed, temperature, and humidity at 10 layers, totaling 30 components. The neural differential equation models the continuous-time dynamical system of the boundary layer, and the neural network includes 3 hidden layers with dimensions of 128, 64, and 32, respectively. The fourth-order Runge-Kutta integrator is used to solve the time evolution process with a time step of 600 seconds. The integrator calculates the forward propagation of the neural network 4 times at each time step, and updates the boundary layer state vector by weighted summation. The solver outputs the turbulent momentum flux and turbulent heat flux distribution at each height layer, as shown in Table 1.

[0156] Table 1 Turbulent flux distribution at different height layers of the boundary layer

[0157]

[0158] Step S09 corrects the future 24-hour wind field prediction results based on the turbulent momentum flux. The technician calculates the wind speed correction amount at each height layer by multiplying the gradient of the turbulent momentum flux with respect to height by the time step of 600 seconds. The average wind speed correction amount in the 10-meter to 50-meter height interval is 0.8 , and the 50-meter to 100-meter interval is 0.6 . The wind speed correction amount is added to the original prediction field to obtain the corrected wind speed vertical distribution. The corrected 10-meter height prediction wind speed is 8.5 50m height is 14.3 The deviations from the actual observations are reduced from 0.7 and 0.9 to 0.2 and 0.3 respectively. The corrected 24h multi-step wind field forecast is used as the final output of the ocean wind field prediction, as shown in Figure 3 .

[0159] The present application brings significant progress compared to the traditional ocean wind field prediction method. The traditional method is based on the dynamic framework of numerical weather prediction models, which needs to solve the atmospheric motion equation set, and has high computational complexity and sensitivity to initial fields. The present application combines Kriging spatial statistics and graph convolutional neural networks to convert the spatial interpolation of sparse observation data into a feature propagation problem on a graph structure, avoiding the high computational cost of traditional physical models. The traditional method has an error accumulation problem in long-term prediction, and the prediction accuracy decreases rapidly as the prediction time increases. The present application uses an encoder-predictor decoupled architecture to extract spatiotemporal invariant feature representations, and through a curriculum learning strategy, the model adapts to its own prediction errors, effectively suppressing the error accumulation phenomenon. The traditional method uses an empirical parameterization scheme for the boundary layer turbulence process, which is difficult to accurately describe the complex turbulence dynamics. The present application models the boundary layer continuous-time dynamical system using neural ordinary differential equations, and learns the turbulence evolution law using a data-driven approach, improving the physical consistency of the vertical distribution prediction of wind speed. In addition, the present application uses a high-order feature interaction modeling algorithm based on tensor decomposition to explicitly capture the nonlinear coupling relationship between multiple observation elements, which can more accurately reconstruct the spatial distribution of the wind field compared to the traditional linear interpolation method. The adaptive rank selection algorithm dynamically adjusts the model complexity according to the observation data density, so that the method maintains stable performance in different sea areas and different observation conditions, and has stronger generalization ability and practical value.

[0160] It should be noted that the variables involved in the present application are explained in detail as shown in Tables 2, 3, and 4.

[0161] Table 2 Variable Explanation Table (First Part)

[0162]

[0163] Table 3 Variable Explanation Table (Second Part)

[0164]

[0165] Table 4 Variable Explanation Table (Third Part)

[0166]

[0167] The above merely illustrates the specific embodiments of the present application, but the protection scope of the present application is not limited thereto, any person skilled in the art can easily think of the changes or replacements within the technical range disclosed by the present application, which should be covered in the protection scope of the present application.

Claims

1. A method for predicting ocean wind fields based on neural networks, characterized in that, This process involves collecting wind speed, wind direction, sea surface temperature, air pressure, and humidity data from marine observation stations to establish a sparse observation dataset. Spatial variogram analysis is performed on this dataset to establish a spatial covariance matrix. This covariance matrix is ​​then converted into a graph-structured adjacency matrix to form a sparse graph topology. The sparse observation dataset and the sparse graph topology are input into a marine wind field spatial reconstruction model. Neighborhood features are aggregated through a graph convolutional layer to output reconstructed wind field data. Historical time series data of the reconstructed wind field are extracted, and spatiotemporal invariant feature representations are extracted using an encoder. These spatiotemporal invariant feature representations are input into a predictor to directly map and generate future multi-step wind field predictions. Wind speed observations are collected to calculate the vertical wind shear rate, determining whether to initiate a boundary layer turbulence correction process. The vertical wind shear rate, sea surface temperature data, and sea surface temperature difference are input into a boundary layer parameterized solver. Turbulent flux is solved using the neural network constant differential equation method. The future multi-step wind field predictions are then corrected based on the turbulent flux.

2. The method according to claim 1, characterized in that, The spatial variability analysis specifically involves calculating the square of the difference between the observations of two ocean observation stations with a distance lag interval, averaging the squares of the observation differences of all ocean observation station pairs that satisfy the distance lag interval, and then dividing by 2 to obtain the semivariance value corresponding to the distance lag interval.

3. The method according to claim 2, characterized in that, The establishment of the spatial covariance matrix specifically involves using the semivariance values ​​corresponding to different distance lags as observation points and fitting a spherical variogram model using the nonlinear least squares method. When the distance is less than the range parameter, the semivariance value of the spherical variogram model increases nonlinearly with the distance, while when the distance is greater than the range parameter, the semivariance value remains at the sill value.

4. The method according to claim 3, characterized in that, The conversion to a graph structure adjacency matrix specifically uses marine observation stations as graph nodes and normalized covariance values ​​as edge weights. When the covariance between two marine observation stations is lower than the covariance threshold, the corresponding edge weight is set to zero.

5. The method according to claim 4, characterized in that, The covariance threshold is determined by calculating the mean and standard deviation of all off-diagonal elements in the spatial covariance matrix, and subtracting 1.5 times the standard deviation from the mean as the covariance threshold.

6. The method according to claim 5, characterized in that, The marine wind field spatial reconstruction model includes an input layer, a first graph convolutional layer, a second graph convolutional layer, a feature interaction layer, an upsampling layer, and an output layer. The first graph convolutional layer performs neighborhood weighted aggregation on the initial node feature vectors according to the graph structure adjacency matrix. The aggregation weight is the product of the edge weight and the attention coefficient.

7. The method according to claim 6, characterized in that, The feature interaction layer adopts a high-order feature interaction modeling algorithm based on tensor decomposition. The node feature tensor composed of the convolutional features of the second graph is decomposed into a product of three factor matrices using CP decomposition. The rank of the CP decomposition is determined by an adaptive rank selection algorithm based on the ocean observation station density parameter, sea area parameter, and historical prediction error parameter.

8. The method according to claim 7, characterized in that, The upsampling layer maps the convolutional features of the second image of the ocean observation station to a regular grid through bilinear interpolation to generate a gridded feature map with a resolution of 0.25 degrees. The output layer contains two convolutional branches that predict the gridded wind speed prediction field and the gridded wind direction prediction field, respectively.

9. The method according to claim 8, characterized in that, The training of the marine wind field spatial reconstruction model uses root mean square error as the loss function for the gridded wind speed prediction field and direction cosine error as the loss function for the gridded wind direction prediction field. The two loss functions are weighted and summed to obtain the total loss function, and gradient descent is performed using the Adam optimizer.

10. The method according to claim 9, characterized in that, The encoder adopts a bidirectional long short-term memory network structure. The forward propagation layer extracts time features from the start time of the historical time series to the future direction, and the backward propagation layer extracts time features from the end time of the historical time series to the past direction. The hidden states in the two directions are concatenated and compressed into a fixed-dimensional spatiotemporal invariant feature representation through a fully connected layer.

Images