Ocean wave multi-element data assimilation method

Abstract

The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-element data assimilation methods, it is related to marine observation data assimilation technical field.The application discloses a kind of ocean sea wave multi-elements data assimilation method, it is related to marine observation data assimilation technical field.

Description

Technical Field

[0001] This invention relates to the field of marine observation data assimilation technology, specifically a method for assimilating multi-element marine wave data. Background Technology

[0002] Multi-element data assimilation of ocean waves is a key method for obtaining high-precision wave reanalysis fields. Existing techniques typically assimilate single elements such as significant wave height, or treat different wave elements as independent variables for separate assimilation, employing methods like optimal interpolation and ensemble Kalman filtering to fuse observational data with numerical model results. These assimilation schemes, when processing multi-element correlation information, only perform simple spatiotemporal matching of each element channel, failing to consider the influence of tidal modulation on the phase state of wave elements at different observation times. Tidal processes cause periodic phase shifts in wave elements across sea surface height and current velocity fields. Observations from different sources do not correspond to consistent tidal phases over time. Directly inputting unaligned multi-element data into the assimilation model introduces spurious inter-element differences, leading to significant phase misalignment and energy distribution distortion in the assimilated wave field nearshore and in strong tidal areas. Furthermore, existing assimilation methods lack in-depth exploration of the coupling relationships between multi-element channels. Complex physical relationships exist among wave elements such as significant wave height, mean wave direction, mean wave period, and sea surface wind speed. For example, the evolution of wave height and wave steepness, wave direction and wind direction, and wave period and wind speed all exhibit nonlinear coupling characteristics. Traditional assimilation models only estimate the error covariance independently for each element channel, making it difficult to capture the energy transfer and mutual constraints across channels. Consequently, the resulting reanalysis field cannot truly reflect the physical consistency among multiple wave elements. To address these issues, it is necessary to solve how to eliminate the phase differences of elements caused by tides during the assimilation process, and how to effectively model and iteratively optimize the inherent coupling between multiple element channels. Summary of the Invention

[0003] This invention provides a method for assimilating multi-element ocean wave data, aiming to eliminate multi-element assimilation errors caused by tidal phase deviations, and to achieve in-depth utilization of the cross-channel coupling relationship of multi-element ocean waves, thereby improving the accuracy of the reanalysis field of ocean wave elements.

[0004] To achieve the above objectives, this invention provides the following technical solution: This invention provides a method for assimilating multi-element ocean wave data. By constructing a multi-element wave state tensor space, multi-source observation data are stacked and expressed on a unified spatiotemporal grid, achieving the fusion and organization of wave elements such as significant wave height, mean wave direction, mean wave period, and sea surface wind speed, providing a standardized multi-dimensional data foundation for subsequent assimilation. A tidal phase modulation function performs phase alignment processing on each element channel. By calculating the phase difference between the instantaneous phase of the element and the tidal phase field and constructing a rotation matrix, rotational projection is performed along the depth and time dimensions, effectively eliminating the phase differences caused by tidal fluctuations among wave elements. This ensures that the dynamic changes among multiple elements remain synchronized in phase, creating favorable conditions for multi-element coupled analysis.

[0005] The multi-factor assimilation model employs three-dimensional convolutional encoding in the encoder, performing block convolution and pooling on the phase-aligned tensor along four dimensions: longitude, latitude, depth, and time, to extract the first latent feature with spatiotemporal compactness. The spectral decomposition layer performs channel-by-channel spectral decomposition on the first latent feature. Preferably, wavelet packet decomposition is used to decompose each feature channel into multiple frequency band subsequences, and based on frequency band energy, the subsequences are distinguished into high-frequency and low-frequency components. This effectively separates rapid disturbances and gradual trends in the wave element, allowing features at different time scales to be analyzed and processed independently. The coupling iteration layer performs cross-channel feature coupling iteration on high-frequency and low-frequency components. It constructs a cross-correlation matrix by calculating the Pearson correlation coefficients of different element channels at the same spatial location and time point, and generates a channel coupling map based on this matrix. Spectral clustering is used to group strongly correlated channels into coupled channel groups, and the features within each group are updated using a weighted average. After multiple iterations, the nonlinear interaction transmission and collaborative correction between elements are achieved, fully exploring the intrinsic physical relationships between wave height, wave direction, period, and wind speed, significantly reducing the impact of single-element observation errors on the overall assimilation results. The decoder performs deconvolution decoding on the updated latent features to generate the assimilated wave element tensor, completing the accurate reconstruction from compressed features to physical space.

[0006] As a technical solution of this invention, when generating the reanalysis field of wave elements for the target sea area based on the assimilated wave element tensor, an inverse transformation corresponding to spectral decomposition is performed on the assimilated data volume of each wave element to restore the spatial distribution field. Furthermore, zero-flux boundary conditions and topography-following boundary conditions are applied based on the coastline boundary and seabed topography to correct boundary constraints, effectively preventing distortion of element values ​​in nearshore and complex topographic regions. During the interpolation to a regular latitude and longitude grid, adaptive neighborhood search and spatial distance-weighted averaging are employed, and the interpolation results are smoothed using Gaussian filtering. This ensures that the generated reanalysis field retains both the small- and medium-scale information of the original observation and model data and possesses spatial continuity and consistency, making it directly applicable to marine forecasting, climate change analysis, and marine engineering.

[0007] The technical effects and advantages provided by the present invention in the above technical solution are as follows: Phase alignment of each element channel in the multi-element wave state tensor space is performed using a tidal phase modulation function. The tidal harmonic constant of the target sea area is obtained to construct the tidal phase field corresponding to the coordinate points in the tensor space. The phase difference between the instantaneous phase of each element channel and the tidal phase field is calculated. Based on the phase difference, a rotation matrix is ​​constructed. This rotation matrix is ​​then used to rotate and project the channel data in the depth and time dimensions, uniformly correcting the data of each element channel to a consistent tidal reference phase. This process eliminates tidal phase differences introduced by different observation times, avoids spurious correlations between elements caused by tidal asynchrony, and ensures that the phase-aligned tensor input to the assimilation model has a unified tidal background. It effectively restores the true spatial and temporal relationships between various wave elements, obtaining a more physically consistent initial assimilation field in strong tidal areas and nearshore waters.

[0008] In the multi-element assimilation model, a spectral decomposition layer and a coupling iteration layer are set up. The spectral decomposition layer performs three-level wavelet packet decomposition on the latent features output by the encoder along each element channel to obtain multiple frequency band subsequences. Based on the frequency band energy, the subsequences are divided into high-frequency subsequences and low-frequency subsequences. After inverse wavelet packet reconstruction, high-frequency components and low-frequency components are formed, realizing the decoupling of multi-element information of ocean waves in the frequency domain. The coupling iteration layer concatenates the high-frequency components and low-frequency components along the element channel dimension, calculates the Pearson correlation coefficient of any two element channels at the same spatial location and time point, constructs a cross-correlation matrix, and establishes a channel coupling graph. By performing spectral clustering based on the normalized Laplacian matrix on the channel coupling graph, closely related element channels are grouped into the same coupling channel group. The channel features within the group are updated by weighted average. The cross-correlation matrix calculation, channel coupling graph construction, spectral clustering, and weighted average update are then iteratively executed. Through multiple iterations, the coupling relationship across channels is gradually transmitted and enhanced, so that the originally isolated ocean wave element features can fully exchange information in the latent space. Without the need to pre-set fixed physical equations, the nonlinear dependence between significant wave height, average wave direction, average wave period, and sea surface wind speed is adaptively captured, so that the assimilated wave element tensor reflects the physical laws of multi-element coordinated change in spatial structure and temporal evolution. The resulting wave element reanalysis field has better physical rationality in dimensions such as joint distribution of wave height and wave direction, and matching of wave period and wind speed. Attached Figure Description

[0009] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0010] Figure 1 This is a flowchart of a method for assimilating multi-element data of ocean waves; Figure 2 This is a flowchart of the construction process for the multi-element wave state tensor space; Figure 3 This is a flowchart of feature processing based on wavelet packet spectral decomposition and cross-channel coupling iteration; Figure 4 This is a schematic diagram of phase alignment of effective wave height time series data. Detailed Implementation

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

[0012] See Figure 1 This invention provides a method for assimilating multi-element ocean wave data, comprising the following steps: constructing a multi-element ocean wave state tensor space, wherein the multi-element ocean wave state tensor space uses longitude, latitude, depth, and time as coordinate axes and multiple ocean wave elements as channels; performing phase alignment processing on each element channel in the multi-element ocean wave state tensor space through a tidal phase modulation function to generate a phase-aligned tensor; inputting the phase-aligned tensor into a multi-element assimilation model, the multi-element assimilation model comprising an encoder, a spectral decomposition layer, a coupled iterative layer, and a decoder, wherein the encoder encodes the phase-aligned tensor to obtain latent features, the spectral decomposition layer performs channel-by-channel spectral decomposition of the latent features to obtain spectral decomposition features, the coupled iterative layer performs cross-channel coupling iteration on the spectral decomposition features to obtain updated latent features, and the decoder decodes the updated latent features to obtain the assimilated ocean wave element tensor; and generating an ocean wave element reanalysis field for the target sea area based on the assimilated ocean wave element tensor.

[0013] In specific implementation, please refer to Figure 2 The method for constructing the multi-element wave state tensor space is as follows.

[0014] Raw wave data from multiple observation sources were acquired for the target sea area. Wave elements include significant wave height, mean wave direction, mean wave period, and sea surface wind speed. Observation sources included satellite remote sensing data, buoy observation data, and numerical model output data. Satellite remote sensing data could be extracted from synthetic aperture radar spectral inversion products; buoy observation data came from real-time records of wave buoys deployed in the target sea area; and numerical model output data was calculated and output by a regional wave numerical model. All acquired raw data were uniformly converted to a common network data format for storage, and the longitude, latitude, depth coordinates, and timestamp of each data point were recorded.

[0015] An unstructured spatial grid is generated by partitioning the target sea area into a spatial grid. The unstructured spatial grid partitioning employs the Delaunay triangulation algorithm, using the target sea area's coastline as the constraint boundary and the seabed topographic slope change rate as the grid density control function. Grid nodes are denser in areas with dramatic topographic changes and sparser in areas with gentle water depth. The spatial grid consists of several unstructured triangular units, with each vertex serving as a grid point. Each grid point has a unique spatial index number and corresponding longitude, latitude, and depth coordinates. The total number of grid points in the unstructured spatial grid is the integer part of the ratio of the target sea area to the preset grid resolution. The preset grid resolution is determined based on the wavelength of the significant wave height spectral peak of the wave numerical model, and is set to one-tenth of the wavelength of the significant wave height spectral peak.

[0016] The observation time range is discretized to establish a unified time series. The observation time range is bounded by the start and end times of interest for the target sea area. The discretized time step is determined based on the shortest revisit period in the satellite remote sensing data, ensuring that at least one type of original data source is updated within each time step. The observation time range is divided into equally spaced time steps, resulting in the time series T = {t1, t2, ..., t...}. M M represents the total number of time points in the time series, and each time point... Each time point corresponds to a timestamp, and the time interval between adjacent time points remains constant.

[0017] The raw data of each wave element are mapped onto a spatial grid and a time series, respectively, to obtain the initial three-dimensional data field for each wave element. For each of the four wave elements—significant wave height, mean wave direction, mean wave period, and sea surface wind speed—the following mapping process is performed.

[0018] All grid points of the spatial grid are selected, and a neighborhood window is constructed centered on each grid point. The shape of the neighborhood window is a cuboid centered on the grid point. The half-width of the window in the longitude direction is the average value of the unstructured grid spatial scale in the longitude direction, and the half-width in the latitude direction is the average value of the unstructured grid spatial scale in the latitude direction. The half-width in the depth direction is 5% of the maximum water depth of the target sea area. No window is set in the time direction at this time because the mapping is performed point-by-point according to the time series. The average values ​​of the spatial scale in the longitude and latitude directions are calculated from the average values ​​of the circumcircle radii of all triangular elements in the spatial grid under the longitude and latitude projections.

[0019] For each time point in the time series, a weighted average of the original data based on the dispersion length is calculated within a neighborhood window to obtain the feature value at each grid point. The calculation method for the weighted average of the dispersion length is as follows: ; in, Indicates at a point in time The feature value at grid point g; K represents the value at the current time point. Below, the number of original data points falling within the neighborhood window of grid point g; v j This represents the feature observation value recorded by the j-th raw data point falling within the neighborhood window. The raw data point originates from satellite remote sensing data, buoy observation data, or numerical model output data; w j Let w represent the weight of the j-th original data point. j Based on the spatial distance d from the j-th original data point to the grid point g j The calculation is performed using the Gaussian kernel function. Spatial distance d j Let be the Euclidean distance between the 3D spatial coordinates of the j-th original data point and the 3D spatial coordinates of the grid point g. The longitude and latitude differences in the Euclidean distance need to be converted to actual geographic distances, while the depth difference is taken as the absolute value. σ is the distance attenuation scale parameter, and σ is the root mean square value of the half-window width in the three spatial directions of the neighborhood window, i.e. , This represents the actual distance corresponding to half the window width in the longitude direction of the neighboring window. This represents the actual distance corresponding to half the window width in the latitudinal direction of the neighboring window. Let K be the half-width of the neighborhood window in the depth direction. If the number of original data points K in the neighborhood window is zero, then double the half-width of the neighborhood window and search again until K is not less than 3. If no original data points are found after expanding the window three times, then... The value is marked as missing and will be filled in by linear interpolation of the valid values ​​from the nearest time points before and after the missing time point in subsequent processing.

[0020] After performing the above-mentioned weighted average of the distances for all grid points, the wave element at the specified time point is obtained. The two-dimensional spatial distribution slices are obtained. The process is repeated for all time points in the time series to obtain the three-dimensional data volume of the wave element. The three dimensions of the three-dimensional data volume are the grid point index, the depth layer index, and the time index, which is the initial three-dimensional data field of the wave element.

[0021] The initial three-dimensional data fields of significant wave height, average wave direction, average wave period, and sea surface wind speed are stacked along the element dimension. This element dimension serves as the fourth dimension, forming a multi-element wave state tensor space. The coordinate axes of this multi-element wave state tensor space are longitude, latitude, depth, and time. The variables corresponding to each coordinate axis are discretized to form a grid point sequence, a depth layer sequence, and a time point sequence. Multiple wave elements are arranged as channels along the channel dimension. Each element in the tensor space uniquely corresponds to a longitude coordinate, latitude coordinate, depth coordinate, time point, and wave element type.

[0022] In practice, the phase alignment of each element channel in the multi-element wave state tensor space is performed by the tidal phase modulation function, and the phase-aligned tensor is generated in the following way.

[0023] Obtain the tidal harmonic constants for the target sea area. The tidal harmonic constants are extracted from the harmonic analysis results output by numerical tidal models or global tidal models. The tidal harmonic constants include the amplitude and lag angle of the major tidal constituents at each grid point in the target sea area. The major tidal constituents include the M2, S2, K1, and O1 tidal constituents. For each grid point in the multi-element wave state tensor space, bilinear interpolation is performed in the tidal harmonic constant dataset based on the longitude and latitude coordinates of the grid point to obtain the amplitude of the M2 tidal constituent at the corresponding grid point. and Chi Jiao The amplitude of the S2 tidal constituents Hechijiao Amplitude of K1 tidal constituents Hechijiao and the amplitude of the O1 tidal fraction. and Chi Jiao Each amplitude and lag angle is stored as a two-dimensional field corresponding to the coordinates of the multi-factor wave state tensor space.

[0024] Based on the tidal harmonic constant, a tidal phase field is constructed corresponding to each coordinate point in the multi-factor wave state tensor space. The construction of the tidal phase field is performed for each spatial coordinate point in the multi-factor wave state tensor space at each time series point. For spatial coordinate point p and time point t... q tidal phase field The calculation method is based on the four main tidal constituents at time point t. q The astronomical tidal phase is a composite phase obtained by amplitude weighting, which is derived by superimposing the complex amplitudes of the tidal constituents and taking their arguments. Tidal phase field. The expression is: ; in, This indicates that spatial coordinate point p is at time point t. q The tidal phase field value at point p is determined by the longitude, latitude, and depth coordinates. This represents the argument function for complex numbers, with the argument taking values ​​in the range [0, 2π). This represents the amplitude of the M2 tidal component at spatial coordinate point p; This represents the amplitude of the S2 tidal component at spatial coordinate point p; This represents the amplitude of the K1 tidal component at spatial coordinate point p; This represents the amplitude of the O1 tidal component at spatial coordinate point p; The angular frequency of the M2 tidal constituent is represented by a value of 1.405189 × 10⁻⁶.-4 radians per second; The angular frequency of the S2 tidal constituent is represented by a value of 1.454441 × 10⁻⁶. -4 radians per second; This represents the angular frequency of the K1 tidal constituent, with a value of 7.292116 × 10⁻⁶. -5 radians per second; The angular frequency of the O1 tidal constituent is represented by a value of 6.759774 × 10⁻⁶. -5 radians per second; t q Representing time point t q The cumulative number of seconds relative to the reference time; This represents the lag angle of the M2 tidal constituent at spatial coordinate point p; This represents the lag angle of the S2 tidal constituent at spatial coordinate point p; The lag angle of the K1 tidal constituent at spatial coordinate point p; The lag angle of the O1 tidal constituent at spatial coordinate point p is represented; e represents the base of the natural logarithm; and i represents the imaginary unit.

[0025] For each feature channel, the phase difference between the instantaneous phase of the feature channel and the tidal phase field is calculated. In the multi-feature wave state tensor space, each feature channel corresponds to the initial three-dimensional data field of a wave feature. For feature channel c, the spatial coordinates p in feature channel c at time point t are extracted. q The feature values ​​at point p are used to construct a time series, which is then subjected to a Hilbert transform to obtain the instantaneous phase of feature channel c at spatial coordinate point p. The Hilbert transform is achieved by performing a Fourier transform on the time series of the element channel c at spatial coordinate point p, multiplying it in the frequency domain by a sign function, and then performing an inverse Fourier transform. Instantaneous phase. The value is taken as the argument of the analytic signal obtained from the Hilbert transform. Phase difference instantaneous phase With tidal phase field The absolute value of the difference, modulo 2π, represents the phase difference. The value of is limited to the interval [0, 2π).

[0026] A rotation matrix is ​​constructed based on the phase difference. For feature channel c at spatial coordinate point p and time point t... q At that point, the rotation matrix R c (p,t q Let R be a 2×2 orthogonal matrix, and let R be the rotation matrix. c (p,t q The elements of ) are determined by phase difference It is composed of the cosine and sine values. Rotation matrix R c (p,t q The first row and first column element of ) has a value of The value of the element in the first row and second column is... The value of the element in the second row and first column is... The value of the element in the second row and second column is... .

[0027] Rotation projection of feature channels along the depth and time dimensions is performed using a rotation matrix to obtain phase-aligned channel data. Rotation projection along the depth and time dimensions means treating the value of feature channel c at spatial coordinate point p along the two-dimensional plane formed by the depth and time dimensions as a two-dimensional vector, and then using the rotation matrix R... c (p,t q This two-dimensional vector undergoes a rotation transformation. In practice, for each fixed longitude and latitude coordinate, the feature channel c can be organized into a two-dimensional matrix along the depth and time dimensions. The rows of the two-dimensional matrix correspond to depth layers, and the columns correspond to time points. For each pair of adjacent depth layers and time points in the two-dimensional matrix, a rotation matrix R is applied. c (p,t q The rotation transformation defined by [reference needed]. The value after rotation and projection replaces the value at the corresponding position in the original feature channel.

[0028] After performing rotational projections of each element channel at all spatial coordinate points and all time points, phase-aligned channel data corresponding to each wave element is obtained. All phase-aligned channel data are then recombine along the channel dimensions. The order of the channel dimensions is consistent with the order of the element channels in the multi-element wave state tensor space, i.e., stacked in the order of significant wave height channel, average wave direction channel, average wave period channel, and sea surface wind speed channel. The resulting four-dimensional data structure is a space composed of longitude, latitude, depth, and time dimensions. The channel dimensions carry the four phase-aligned wave element data, generating a phase-aligned tensor.

[0029] See Figure 4The horizontal axis of the graph represents time in hours, covering a 50-hour observation period. The vertical axis represents the phase quantity of the significant wave height sequence after normalization and tidal phase field processing, in units of significant wave height (meters) divided by the normalized tidal phase. The blue solid line represents the original significant wave height time series without phase alignment processing. Its curve shows a relatively flat change with small amplitude over 48 hours, with a maximum value of about 3.5 meters. The fluctuation period is not obvious, and the peaks and troughs are relatively evenly distributed. The red dashed line represents the significant wave height sequence after phase alignment processing using the tidal phase modulation function. Compared to the original significant wave height sequence, the red dashed line curve shows obvious periodic oscillation characteristics, periodically corresponding to the changes in tidal phase. The amplitude ranges from about -4 meters to 4 meters, with alternating peaks and troughs and a significantly increased amplitude, indicating that the phase alignment processing effectively enhances the tidal correlation of the significant wave height sequence. The green dotted line represents the normalized tidal phase field, which is relatively stable, basically maintaining in the range of 2 to 2.5, showing a slow sinusoidal fluctuation trend. The overall curve change process can be divided into multiple tidal cycle stages, with a complete cycle of peaks and troughs occurring approximately every 6 hours. This figure fully demonstrates the effect of using the tidal phase modulation function to perform phase alignment processing on the effective wave height channel in the multi-element wave state tensor space in Example 2. Through phase alignment, the effective wave height sequence is synchronized with the tidal phase, improving the spatiotemporal coherence of the multi-element wave data and providing a phase-consistent wave data foundation for subsequent assimilation model input.

[0030] In practice, the phase-aligned tensor is input into the multi-factor assimilation model, which includes an encoder, a spectral decomposition layer, a coupled iterative layer, and a decoder, with each module connected sequentially.

[0031] The input to the multi-factor assimilation model is a phase alignment tensor, which has longitude, latitude, depth, time, and channel dimensions. The channel dimension corresponds to four wave elements. The encoder receives the phase alignment tensor, performs 3D convolutional encoding on it, and outputs the first latent feature. The specific method of the encoder performing 3D convolutional encoding is as follows.

[0032] The phase alignment tensor is divided into blocks along four dimensions: longitude, latitude, depth, and time, resulting in multiple spatial-temporal blocks. During block division, the block size along the longitude dimension is 10 grid points, with a step size of 5 grid points between blocks; the block size along the latitude dimension is 10 grid points, with a step size of 5 grid points between blocks; the block size along the depth dimension covers all depth layers and is not divided; the block size along the time dimension is 4 time points, with a step size of 2 time points between blocks. The generated spatial-temporal blocks overlap in the longitude and latitude dimensions, with an overlap ratio of 50% of the block size. The dimensions of each spatial-temporal block are: number of longitude grid points within the block × number of latitude grid points within the block × number of depth layers × number of time points within the block × number of channels.

[0033] Each spatial-temporal block is input into the first 3D convolutional layer. This first 3D convolutional layer contains eight convolutional kernels, each with dimensions of 3, 3, 3, and 2 along the longitude, latitude, depth, and time dimensions, respectively. The stride of each kernel in all four dimensions is set to 1, and zero values ​​are padded at the edges of each of the four dimensions (longitude, latitude, depth, and time) to ensure the output dimensions match the size of the input spatial-temporal block. The first 3D convolutional layer operates independently on each spatial-temporal block, with each kernel performing a four-dimensional convolution on the block, outputting a feature map. The eight kernels output a total of eight feature maps, each with dimensions identical to the four spatial-temporal dimensions of the block. The kernel weights of the first 3D convolutional layer are initialized using a truncated normal distribution with a mean of 0 and a standard deviation of 0.02, and the bias is initialized to 0.

[0034] The output of the first 3D convolutional layer is subjected to batch normalization and non-linear activation. The batch normalization operation calculates the mean and variance in both the batch dimension and the spatial-temporal dimension, normalizing each feature map. The momentum parameter in batch normalization is set to 0.99, and the minimum constant to prevent division by zero is set to 1×10⁻⁻⁶. 5 After batch normalization, a nonlinear activation function is applied. This nonlinear activation function uses a leaky modified linear unit (LMU), with a slope of 0.2 in the negative region. The LMU independently calculates the activation value for each element in each feature map. After batch normalization and nonlinear activation, the first convolutional feature is obtained, containing eight feature maps. The spatial-temporal dimensions of these feature maps are the same as those of the spatial-temporal blocks.

[0035] The first convolutional feature is input into the second 3D convolutional layer. The second 3D convolutional layer contains 16 convolutional kernels, each with dimensions of 3, 3, 1, and 2 along the longitude, latitude, depth, and time dimensions, respectively. The kernels in the second 3D convolutional layer have a stride of 1 in all four dimensions, and zero values ​​are padded at the edges of the longitude, latitude, depth, and time dimensions to maintain a constant output size. The second 3D convolutional layer receives the first convolutional feature, performs a four-dimensional convolution on its eight feature maps, and outputs 16 feature maps. The kernel weights of the second 3D convolutional layer are initialized using the same method as those of the first 3D convolutional layer.

[0036] Max pooling is performed on the output of the second 3D convolutional layer. The pooling window for max pooling is 2, 2, 1, and 2 along the longitude, latitude, depth, and time dimensions, respectively. The stride of the pooling window is 2 in the longitude direction, 2 in the latitude direction, 1 in the depth direction, and 2 in the time direction. Max pooling takes the maximum value within each pooling window, reducing the spatial and temporal dimensions of each feature map output by the second 3D convolutional layer while keeping the depth dimension unchanged. The features output by max pooling are used as the first latent features. The first latent features consist of 16 feature maps, with dimensions of the reduced longitude, latitude, depth, and time. The feature channel dimension of the first latent features is composed of 16 feature maps, each of which is considered a feature channel.

[0037] The spectral decomposition layer receives the first latent feature. It performs channel-by-channel spectral decomposition on each element channel of the first latent feature, decomposing each element channel into high-frequency and low-frequency components. Channel-by-channel spectral decomposition performs multi-level wavelet packet decomposition on the time-dimensional data of each element channel, dividing the time-dimensional signal into multiple frequency bands. Based on the frequency band energy, each band is further divided into high-frequency and low-frequency components, which are then reconstructed to obtain high-frequency and low-frequency tensors respectively. The spectral decomposition layer stacks the high-frequency tensors of all element channels along the element channel dimension to obtain the high-frequency components; it also stacks the low-frequency tensors of all element channels along the element channel dimension to obtain the low-frequency components. The number of element channels for both the high-frequency and low-frequency components is the same as the number of element channels for the first latent feature.

[0038] The coupling iteration layer receives high-frequency and low-frequency components. It performs cross-channel feature coupling iteration on these components. The layer concatenates the high-frequency and low-frequency components along the feature channel dimension to obtain the feature to be coupled. It calculates the Pearson correlation coefficient between feature values ​​at the same spatial location and time point for any two feature channels in the feature to be coupled, forming a cross-correlation matrix. Based on the cross-correlation matrix, it constructs a channel coupling graph and performs spectral clustering on the graph to obtain coupled channel groups. It then performs a weighted average of the channel features within each coupled channel group and updates the features. This process is iteratively executed until the preset number of iterations is reached. After iteration, the coupling iteration layer reassembles the iterated high-frequency and low-frequency components along the feature channel dimension to obtain the second latent feature. The number of feature channels in the second latent feature is the same as that in the first latent feature.

[0039] The decoder receives the second latent feature. It performs deconvolution decoding on the second latent feature to obtain the assimilated wave element tensor. The decoder consists of two cascaded deconvolutional layers and an output convolutional layer. The first deconvolutional layer contains 16 deconvolutional kernels, each with dimensions of 3, 3, 1, and 2 along the longitude, latitude, depth, and time dimensions, respectively. The stride during deconvolution is set to 2, 2, 1, and 2. The deconvolution operation doubles the size of the second latent feature in the longitude, latitude, and time dimensions, while keeping the depth dimension unchanged. The second deconvolutional layer contains 8 deconvolutional kernels, each with dimensions of 3, 3, 3, and 2 along the longitude, latitude, depth, and time dimensions, respectively. The stride is set to 1, 1, 1, and 1. The deconvolution operation restores the original size of the spatial-temporal block. The output convolutional layer contains four convolutional kernels, each with a size of 1 along the longitude, latitude, depth, and time axes, and a stride of 1. The output convolutional layer maps the eight feature maps from the second deconvolutional layer into four channels, corresponding to significant wave height, average wave direction, average wave period, and sea surface wind speed. The spatial-temporal dimension of the assimilated wave element tensor output by the decoder is consistent with the spatial-temporal dimension of the input block. The assimilated wave element tensors corresponding to all spatial-temporal blocks are then weighted and concatenated according to the overlapping regions during block segmentation. The weights of the overlapping regions are applied using a linear gradient-in / gradual-out method, ultimately yielding an assimilated wave element tensor covering the entire target sea area.

[0040] In practice, the multi-factor assimilation model needs to be trained before use. The training process is as follows.

[0041] A training dataset is constructed. The input samples of the training dataset are phase-aligned tensors obtained by processing the historical multi-factor wave state tensor space using a tidal phase modulation function. The labels are high-resolution wave reanalysis products from the same spatiotemporal range that have undergone quality control and manual verification. The training samples are organized in a spatial-temporal block manner, with the spatial-temporal block size consistent with the encoder block size. The training dataset includes a training set, a validation set, and a test set, with a partitioning ratio of 8:1:1.

[0042] During the training of the multi-factor assimilation model, an adaptive moment estimation optimization algorithm was used. The initial learning rate of the adaptive moment estimation optimization algorithm was set to 0.001, and the first-order moment decay factor was set to... Set to 0.9, second-order moment decay factor Set to 0.999, numerical stability constant Set to 10 -8 The training batch size is set to 16, and one training epoch consists of iterating through all training samples. The total number of training epochs is set to 200. The loss function is weighted mean squared error, and its expression is as follows: ; Where L represents the loss value during the training process of the multi-factor assimilation model; This represents the number of spatial-temporal block samples in a training batch. This indicates the number of wave element channels, with a value of 4. This represents the loss weighting coefficient for the c-th wave element channel. Determined based on the observation accuracy and rate of change of physical quantities of each wave element; for the significant wave height channel, The value is 1.0; for the average waveguide channel, The value is 0.5; for the average wave period channel, The value is 0.8; for the sea surface wind speed channel, The value is set to 0.7; the above weighting coefficient setting ensures that the contribution of each wave element in the loss function is on the same order of magnitude, avoiding the dominance of gradient updates by elements with large values ​​such as effective wave height and sea surface wind speed. This represents the label value of the r-th spatial-temporal block sample in the c-th wave element channel in the training batch. The spatiotemporal range of the label value corresponds completely to the spatial-temporal block sample. This represents the assimilation value output by the multi-factor assimilation model for the r-th spatial-temporal block sample in the c-th wave element channel; r represents the index number of the spatial-temporal block sample within the training batch; and c represents the index number of the wave element channel.

[0043] During training, the backpropagation algorithm is used to calculate the gradient of the loss function L with respect to all learnable parameters in the multi-factor assimilation model. These learnable parameters include the kernel weights and biases of the first 3D convolutional layer, the kernel weights and biases of the second 3D convolutional layer, the deconvolution kernel weights and biases of the first and second deconvolutional layers, the kernel weights and biases of the output convolutional layer, the scaling and offset parameters in batch normalization, and the weight parameters involved in the coupled iterative layers. Wavelet packet decomposition in the spectral decomposition layer uses fixed filter coefficients and is not involved in training. Every 5 training epochs, the validation loss is calculated on the validation set, and its changes are monitored. If the validation loss does not reach a new low for 20 consecutive training epochs, the learning rate is reduced to 0.5 times the current learning rate. The training continues until the learning rate drops below 0.001 times the initial learning rate or the validation loss fluctuates by less than 10% over 30 consecutive training epochs. -6 Training terminates at the specified time, and the model parameters corresponding to the training epoch with the lowest validation loss are selected as the final parameters of the multi-factor assimilation model. The trained multi-factor assimilation model is then used for the assimilation of wave elements in the target sea area.

[0044] In specific implementation, please refer to Figure 3 The process of performing channel-by-channel spectral decomposition on the first latent feature to obtain high-frequency and low-frequency components is as follows.

[0045] The spectral decomposition layer receives the first latent feature output from the encoder. This first latent feature contains 16 feature channels, each with a data structure that is a multidimensional array along four dimensions: longitude, latitude, depth, and time. Channel-by-channel spectral decomposition is performed independently for each feature channel of the first latent feature, with the processing flow being identical for each channel. Taking the u-th feature channel as an example, the channel-by-channel spectral decomposition process is as follows: extract the time series of all spatial location points in the u-th feature channel. These spatial location points are determined by the longitude grid point index, latitude grid point index, and depth layer index. Each spatial location point corresponds to a time series of length M. t Time series, M t This represents the total number of time points in the time series. For each time series, a three-level wavelet packet decomposition is performed along the time axis.

[0046] The three-level wavelet packet decomposition uses the db4 wavelet basis function from the Daubechies wavelet family. The db4 wavelet basis function has fourth-order vanishing moments and a filter length of 8. Each level of wavelet packet decomposition passes the input signal through a pair of orthogonal mirror filters. The orthogonal mirror filters consist of a low-pass decomposition filter and a high-pass decomposition filter. The coefficient sequence of the low-pass decomposition filter is composed of the scaling function coefficients of the db4 wavelet basis function, and the coefficient sequence of the high-pass decomposition filter is composed of the wavelet function coefficients of the db4 wavelet basis function. In the first-level wavelet packet decomposition, the input time series is convolved with the low-pass and high-pass decomposition filters respectively, and the convolution results are downsampled at intervals with a downsampling factor of 2 to obtain the first-level low-frequency subsequence and the first-level high-frequency subsequence. In the second-level wavelet packet decomposition, the first-level low-frequency subsequence is passed through a low-pass decomposition filter and a high-pass decomposition filter, and then downsampled at intervals to obtain the first part and the second part of the second-level low-frequency subsequence. Similarly, the first-level high-frequency subsequence is passed through a low-pass decomposition filter and a high-pass decomposition filter, and then downsampled at intervals to obtain the first part and the second part of the second-level high-frequency subsequence. In the third-level wavelet packet decomposition, the four subsequences obtained from the second-level decomposition are again passed through a low-pass decomposition filter and a high-pass decomposition filter, and then downsampled at intervals, ultimately resulting in eight frequency band subsequences. These eight frequency band subsequences are arranged sequentially from low to high frequency on the frequency axis, with the first frequency band subsequence occupying the entire frequency range. Partially, the 8th frequency band subsequence occupies the entire frequency range. part.

[0047] For each frequency band subsequence, calculate the frequency band energy of the frequency band subsequence. The calculation method is the sum of squares of each sampling point in the subsequence of that frequency band: ; in, This represents the frequency band energy of the b-th frequency band subsequence of the u-th feature channel at the s-th spatial location point; u represents the index number of the feature channel, with a value ranging from 1 to 16; s represents the index number of the spatial location point, which is determined by a linear combination of the longitude grid point index, the latitude grid point index, and the depth layer index, with a value ranging from 1 to 16. , is the total number of all spatial location points; b represents the index number of the frequency band subsequence. After three-level wavelet packet decomposition, a total of 8 frequency band subsequences are generated. The value of b ranges from 1 to 8. b=1 corresponds to the lowest frequency frequency band subsequence, and b=8 corresponds to the highest frequency frequency band subsequence. This represents the number of sampling points in the b-th frequency band subsequence. Since each level of decomposition involves downsampling at intervals, the number of sampling points in each frequency band subsequence after the three-level wavelet packet decomposition is equal to the length M of the original time series. t Round down to one-eighth of the total. This represents the wavelet packet coefficient value of the u-th element channel at the b-th frequency band subsequence at the s-th spatial location point; q represents the index number of the sampling point in the frequency band subsequence.

[0048] Calculate the average band energy of the u-th element channel at the s-th spatial location point. The average band energy is the sum of the band energies of the 8 band sub-sequences at that spatial location point divided by 8. Label the band sub-sequences whose band energy at the s-th spatial location point exceeds the average band energy as high-frequency sub-sequences, and label the band sub-sequences whose band energy does not exceed the average band energy as low-frequency sub-sequences.

[0049] Inverse wavelet packet reconstruction is performed on all frequency band subsequences marked as high-frequency subsequences. This inverse wavelet packet reconstruction is the inverse operation of wavelet packet decomposition. Wavelet packet coefficients in the frequency band subsequences marked as high-frequency subsequences are retained, while all wavelet packet coefficients in the frequency band subsequences marked as low-frequency subsequences are set to zero. For the eight frequency band subsequences after zeroing, upsampling and filtering reconstruction are performed progressively from the third level. Upsampling is achieved by inserting zero values ​​between each sampling point, with an upsampling factor of 2. Filtering reconstruction uses the reconstruction low-pass filter and reconstruction high-pass filter corresponding to the db4 wavelet basis function. The coefficient sequence of the reconstruction low-pass filter is the reverse order of the coefficient sequence of the decomposed low-pass filter, and the coefficient sequence of the reconstruction high-pass filter is the reverse order of the coefficient sequence of the decomposed high-pass filter. During the reconstruction process, each pair of low-frequency and high-frequency subsequences is passed through the reconstruction low-pass filter and reconstruction high-pass filter respectively, and then summed to progressively restore the original time series length. The reconstructed sequence is the high-frequency reconstructed time series of the u-th element channel at the s-th spatial location point. Arrange the high-frequency reconstructed time series of all spatial locations according to their original spatial locations to obtain the high-frequency tensor corresponding to the u-th feature channel. The dimension of the high-frequency tensor is consistent with the dimension of the u-th feature channel in the first latent feature.

[0050] Inverse wavelet packet reconstruction is performed on all frequency band subsequences labeled as low-frequency subsequences. The reconstruction method is the same as that for high-frequency reconstruction, except that all wavelet packet coefficients in the frequency band subsequences labeled as high-frequency subsequences are set to zero, while the wavelet packet coefficients in the frequency band subsequences labeled as low-frequency subsequences are retained. The reconstructed sequence is the low-frequency reconstructed time series of the u-th feature channel at the s-th spatial location point. The low-frequency reconstructed time series of all spatial location points are arranged according to their original spatial locations to obtain the low-frequency tensor corresponding to the u-th feature channel.

[0051] The above operation is performed on all 16 feature channels of the first latent feature, resulting in 16 high-frequency tensors and 16 low-frequency tensors. The 16 high-frequency tensors are concatenated along the feature channel dimension; the concatenated result is the high-frequency component, which has 16 feature channels. Similarly, the 16 low-frequency tensors are concatenated along the feature channel dimension; the concatenated result is the low-frequency component, which also has 16 feature channels. The high-frequency and low-frequency components are identical in size across the spatial, temporal, and feature channel dimensions.

[0052] The coupling iteration layer receives high-frequency and low-frequency components and performs cross-channel feature coupling iteration. The coupling iteration process involves concatenating the high-frequency and low-frequency components along the feature channel dimension. The high-frequency component has 16 feature channels, and the low-frequency component also has 16 feature channels, resulting in a total of 32 feature channels for the feature to be coupled after concatenation. For any two feature channels in the feature to be coupled, the Pearson correlation coefficient is calculated between feature values ​​at the same spatial location and time point. For feature channel a and feature channel b, the Pearson correlation coefficient is calculated. The calculation method is as follows: Flatten the feature values ​​of all spatial locations and all time points in feature channel a into a one-dimensional vector, and flatten the feature values ​​of all spatial locations and all time points in feature channel b into a one-dimensional vector. Calculate the Pearson correlation coefficient between the two one-dimensional vectors. Repeat this calculation for all possible feature channel pairs to obtain a 32-row, 32-column cross-correlation matrix. The element in the a-th row and b-th column of the cross-correlation matrix is ​​the Pearson correlation coefficient between feature channel a and feature channel b. The cross-correlation matrix is ​​a symmetric matrix, and the diagonal elements are always 1.

[0053] A channel coupling graph is constructed based on the cross-correlation matrix. The channel coupling graph has a total of 32 nodes, each corresponding to one feature channel to be coupled. The edge weight between node a and node b in the channel coupling graph is the absolute value of the Pearson correlation coefficient between the corresponding feature channels. Set edge weight thresholds. ,like In the channel coupling graph, the edge between node a and node b is not retained, and the edge weight threshold is set. The basis for setting it to 0.3 is that a Pearson correlation coefficient below 0.3 indicates that there is only a weak correlation between the two variables, which contributes little to cross-channel coupling. Discarding it can reduce noise interference to channel grouping.

[0054] Perform spectral clustering based on the normalized Laplacian matrix on the channel coupling graph. During spectral clustering, construct the adjacency matrix W of the channel coupling graph, and the elements of the adjacency matrix W... when ,otherwise The degree matrix D is calculated from the adjacency matrix W. The degree matrix D is a diagonal matrix, and the diagonal elements are... Calculate the normalized Laplace matrix. For the normalized Laplace matrix Perform eigenvalue decomposition and take the minimum value. The eigenvectors corresponding to each eigenvalue. This is the preset number of clusters. The value is determined based on the number of connected components in the channel coupling graph. These connected components are found in the subgraph that retains connections above a certain edge weight threshold using a depth-first search algorithm. The number of connected components is then used as the threshold value. The value of . All 32 nodes corresponding to The feature matrix is ​​composed of eigenvectors. K-means clustering is performed on each row of the feature matrix, with the number of centers in each k-means cluster set to 1. The clustering results divided the 32 element channels into One coupling channel group.

[0055] A weighted average is performed on the channel features within each coupled channel group. For a coupled channel group containing η feature channels, the feature values ​​of all feature channels within the coupled channel group at the same spatial location and time point are weighted and averaged. The weight of the average is the average edge weight between each feature channel and other feature channels within the coupled channel group. The value of the weighted average replaces the original feature values ​​of all feature channels within the coupled channel group at the corresponding spatial location and time point. After performing the weighted average update on all coupled channel groups, the updated channel features are obtained.

[0056] The updated channel features are then used again as the features to be coupled, and the cross-correlation matrix calculation, channel coupling graph construction, spectral clustering, and weighted average update operations are performed iteratively. The preset number of iterations is set to 5, which is determined based on the convergence characteristics of the cross-correlation matrix during the coupling iteration process. After 5 iterations, the change in the Frobenius norm of the cross-correlation matrix usually decreases to less than 5% of the initial change. In the 5 iterations, the input for the first iteration is the feature to be coupled obtained by concatenating the high-frequency and low-frequency components; the input for the second to fifth iterations is the channel features updated after the previous iteration.

[0057] After the 5th iteration, no further iterations are performed. Following the concatenation order of high-frequency and low-frequency components before the 1st iteration, the channel features updated in the 5th iteration are split along the element channel dimension. During splitting, the updated features of the first 16 element channels are taken as the updated high-frequency components, and the updated features of the last 16 element channels are taken as the updated low-frequency components. The updated high-frequency and low-frequency components are then recombine along the element channel dimension. This recombination is achieved by concatenating the updated high-frequency and low-frequency components along the element channel dimension to generate a second latent feature. The second latent feature has 32 element channels, and its spatial and temporal dimensions are the same as the first latent feature.

[0058] In practice, the process of generating the reanalysis field of wave elements for the target sea area based on the assimilated wave element tensor is as follows.

[0059] The assimilated wave element tensor is split along the channel dimension. The channel dimension of the assimilated wave element tensor contains four wave elements, arranged in the channel order agreed upon during the encoding phase: significant wave height channel, mean wave direction channel, mean wave period channel, and sea surface wind speed channel. The splitting operation along the channel dimension involves dividing the assimilated wave element tensor by index along the channel dimension, resulting in four independent data volumes. Each independent data volume is the assimilated data volume for one wave element. The dimensions of the assimilated data volume are longitude, latitude, depth, and time. The assimilated data volume stores the assimilated values ​​of the corresponding wave element on the target sea area's three-dimensional spatial grid and a unified time series.

[0060] For each wave feature's assimilated data volume, an inverse transform corresponding to the channel-by-channel spectral decomposition performed by the spectral decomposition layer is executed. In the spectral decomposition layer, channel-by-channel spectral decomposition performs a third-level wavelet packet decomposition on each feature channel of the first latent feature, obtaining high-frequency and low-frequency components. The corresponding inverse transform is a third-level wavelet packet reconstruction, which performs a complete third-level wavelet packet inverse transform on the time series of each spatial location point in the assimilated data volume. For the assimilated data volume of significant wave height, the processing method involves extracting the time series at each spatial location point in the significant wave height assimilated data volume. The spatial location points are determined by the longitude grid point index, latitude grid point index, and depth layer index. A third-level wavelet packet decomposition is performed on the extracted time series, using the same wavelet basis function used in the spectral decomposition layer, namely the db4 wavelet basis function, with the filter coefficients remaining unchanged. The third-level wavelet packet decomposition decomposes the time series into eight frequency band subsequences. Unlike the thresholding process in the spectral decomposition layer, the inverse transform retains the wavelet packet coefficients of all eight frequency band subsequences without performing any coefficient zeroing operations. Three-level wavelet packet reconstruction was performed on all eight frequency band subsequences. The reconstruction filters were the reconstruction filter bank corresponding to the db4 wavelet basis functions. The reconstruction process involved upsampling and filtering synthesis in sequence. After reconstruction, the original time series length was restored, resulting in the reconstructed time series at that spatial location. The reconstructed time series completely preserved all frequency components of the original signal in the frequency domain. However, after feature updates through the spectral decomposition layer and the coupled iterative layer, the reconstructed values ​​reflected the state after assimilation of multiple factors.

[0061] The inverse transformation described above is performed on all spatial locations, and the reconstructed time series is arranged according to spatial location. The reconstructed time series forms the spatial distribution field of this wave element. The spatial coordinate framework of the spatial distribution field is consistent with the assimilated data volume, that is, it adopts a three-dimensional coordinate system of unstructured spatial grid. The spatial distribution fields of significant wave height, mean wave direction, mean wave period, and sea surface wind speed are all generated in the above manner.

[0062] Based on the coastline boundaries and seabed topography data of the target sea area, the spatial distribution field of each wave element is corrected by boundary constraints to obtain the corrected spatial distribution field. Boundary constraint correction includes two aspects of constraint processing.

[0063] The method for correcting coastline boundary constraints is as follows: Obtain vector data of the coastline boundary of the target sea area. This vector data includes the longitude and latitude coordinates of the coastline inflection point sequence. Project the coastline boundary onto the longitude-latitude plane of the spatial distribution field and identify the grid points located on the coastline boundary. The method for determining grid points on the coastline boundary is to calculate the shortest distance between the horizontal projection point of each grid point and the broken line segment of the coastline. If the shortest distance is less than half the average side length of the spatial grid, the grid point is marked as a coastline boundary grid point. Set the feature values ​​at the coastline boundary grid points as zero-flux boundary conditions. The mathematical expression of the zero-flux boundary condition is that the feature gradient in the direction of the coastline normal is zero. On the discrete grid, the zero-flux boundary condition is implemented by assigning the feature values ​​of the coastline boundary grid points the feature values ​​of one grid point inward along the normal direction.

[0064] The seabed topography boundary constraint correction method is as follows: Seabed topography data for the target sea area is acquired. This data is regular gridded water depth data, with depth values ​​in meters. For each grid point in the spatial distribution field, the depth coordinates of the grid point are compared with the seabed topography data. The depth coordinates use the sea surface as the zero level, with negative values ​​downwards. When the depth coordinate value of a grid point is less than or equal to the negative absolute value of the water depth at the corresponding horizontal position in the seabed topography data, the grid point is located at or below the seabed topography surface. Grid points located at or below the seabed topography surface are marked as seabed topography boundary grid points. Feature values ​​at the seabed topography boundary grid points are set to a topography-following boundary condition. The topography-following boundary condition allows feature values ​​to exist along the seabed topography surface, but the gradient of the feature values ​​along the topography normal direction is zero. On the discrete grid, the topography-following boundary condition sets the feature values ​​of the seabed topography boundary grid points to the feature values ​​of the grid points one grid point above the seabed normal direction. If no grid point is available in the upward direction, the feature value is set to zero. In areas where the seabed topographic slope exceeds 30 degrees, the topographic follow-boundary condition is degraded to a no-slip boundary condition, that is, the feature value is directly set to zero, in order to avoid numerical anomalies caused by interpolation in steep areas.

[0065] The corrected spatial distribution fields of all wave elements are combined. The combination operation involves stacking the corrected spatial distribution fields of significant wave height, mean wave direction, mean wave period, and sea surface wind speed along the element dimensions in the same channel order to generate a multi-element corrected data volume. The spatial coordinate system of the multi-element corrected data volume is still a four-dimensional data structure under an unstructured spatial grid.

[0066] The multi-element corrected data volume is interpolated onto a regular latitude and longitude grid of the target sea area to generate a reanalysis field of ocean wave elements. The regular latitude and longitude grid is established by setting the longitude resolution of the grid based on the longitude and latitude span of the target sea area. and latitude resolution Longitude resolution and latitude resolution All resolutions are set to 0.1 degrees. This 0.1-degree setting is based on the spatial characteristic scale of mesoscale eddies in the target sea area. A 0.1-degree grid resolution can distinguish wave wave characteristics with wavelengths greater than 20 kilometers. The grid points of the regular latitude and longitude grid are located at the intersections of meridians and parallels of latitude. The total number of grid points in the regular latitude and longitude grid is determined by the longitude and latitude ranges of the target sea area and the resolution.

[0067] The interpolation operation onto the regular latitude and longitude grid is performed separately for the corrected spatial distribution field of each wave element. Taking the corrected spatial distribution field of significant wave height as an example, the interpolation process is as follows: obtain the discrete point coordinates and discrete point values ​​of the corrected spatial distribution field of significant wave height. The discrete point coordinates are the three-dimensional coordinates of the grid points in the corrected spatial distribution field of significant wave height in the unstructured spatial grid, including longitude, latitude, and depth values. The discrete point values ​​are the corrected significant wave height values ​​at the corresponding grid points.

[0068] Using each grid point of the regular latitude and longitude grid as the center, discrete points within a horizontal search radius of 50 kilometers are used as the interpolation neighborhood. The preset search radius of 50 kilometers is determined based on the following: the horizontal coherence scale of mesoscale ocean processes in the target sea area is typically between 30 and 80 kilometers, and a 50-kilometer search radius can cover the variability characteristics of wave elements at a medium spatial scale. During the search, the longitude and latitude values ​​of the regular latitude and longitude grid points and each discrete point are converted into spherical distances. The spherical distance is calculated using the great circle distance formula, specifying the Earth's average radius as 6371 kilometers. For longitude differences... and latitude difference The spherical distance S is calculated as follows: ; Where S represents the spherical distance from a regular latitude and longitude grid point to a discrete point, in kilometers; 6371 represents the average radius of the Earth, in kilometers; and arccos represents the inverse cosine function. The latitude value of the regular latitude and longitude grid points is represented in radians. Represents the latitude value of a discrete point, in radians; The longitude value of the regular latitude and longitude grid points is represented in radians. Represents the longitude value of a discrete point, in radians.

[0069] When the spherical distance S is less than or equal to 50 kilometers, the corresponding discrete point is selected into the interpolation neighborhood.

[0070] For discrete points within the interpolation neighborhood, the spatial distance from each discrete point to a grid point is used as the weight, with smaller spatial distances resulting in larger weights. A weighted average of the discrete point values ​​is calculated to obtain the initial interpolation at each grid point. The weighting function for the weighted average uses the inverse squared distance form, where the weight of the d-th discrete point is... Defined as: ; in, This represents the weight of the d-th discrete point within the interpolation neighborhood; This represents the spherical distance from the d-th discrete point in the interpolation neighborhood to a point on the regular latitude and longitude grid, in kilometers. To prevent the introduction of a minimal constant when the denominator is zero, The value is 10 -6 Square kilometers.

[0071] Initial interpolation at grid points The weighted average of the values ​​of all discrete points within the interpolation neighborhood: ; in, This represents the initial interpolation at point g in the regular latitude and longitude grid. This represents the total number of discrete points within the interpolation neighborhood corresponding to the regular latitude and longitude grid point g; This represents the weight of the d-th discrete point within the interpolation neighborhood; This represents the discrete point value of the d-th discrete point within the interpolation neighborhood; d represents the index number of the discrete point within the interpolation neighborhood.

[0072] If the number of discrete points in the interpolation neighborhood of regular latitude and longitude grid points If there are fewer than 3, increase the search radius to 75 kilometers and search again. If the radius is increased... If the number is still less than 3, the search radius will be further increased to 100 kilometers. If the search radius is 100 kilometers... If there are still fewer than 3, the regular latitude and longitude grid points will be marked as land or data blank areas, and the initial interpolation will be set to null values, which will not participate in the smoothing calculation in subsequent processing.

[0073] The initial interpolation is smoothed using a Gaussian filter to obtain smoothed grid point values. The Gaussian filtering smoothing is performed on a two-dimensional plane formed by the longitude and latitude directions of a regular latitude and longitude grid. The kernel size of the Gaussian filter is 3×3 along both the longitude and latitude directions, meaning the Gaussian filter kernel covers the grid point itself and its eight surrounding adjacent grid points. The weight matrix of the Gaussian filter kernel is obtained by sampling a two-dimensional Gaussian function onto discrete grid points. The standard deviation of the two-dimensional Gaussian function is also included. The value is set to 1 grid spacing, with a longitude resolution of 0.1 degrees and a latitude resolution of 0.1 degrees. The corresponding spatial scale is approximately 10 kilometers. Gaussian filtering smoothing removes short-wavelength oscillations that may be caused by interpolation from unstructured to regular grids. The smoothed grid point values ​​are shown below. The calculation method is to perform discrete convolution with a Gaussian filter kernel and initial interpolation within a 3×3 window centered on grid point g. During convolution, empty grid points within the window are skipped, and empty grid points are not included in the convolution or assigned weights. At the same time, the weights of non-empty grid points are normalized.

[0074] The smoothed grid point values ​​are used as the reanalysis field values ​​of the significant wave height on the regular latitude and longitude grid. The same interpolation process is performed on the corrected spatial distribution fields of the mean wave direction, mean wave period, and sea surface wind speed, respectively, to obtain the reanalysis field values ​​of the wave elements for the mean wave direction, mean wave period, and sea surface wind speed. During the interpolation of the mean wave direction, since the mean wave direction is an angular quantity, the angle values ​​need to be converted into unit direction vectors during weighted averaging. The horizontal and vertical components of the unit direction vector are then weighted and averaged separately, and the synthesized vector direction is converted back into angle values. These converted angle values ​​are used as the values ​​of the mean wave direction on the regular latitude and longitude grid. After the corrected spatial distribution field of all wave elements is interpolated and correlated with the depth coordinates, a multi-element wave element reanalysis field for the target sea area is generated. The data format of the wave element reanalysis field is a common network data format, which includes longitude coordinate variables, latitude coordinate variables, time coordinate variables, and reanalysis numerical variables of each wave element.

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

Claims

1. A method for assimilating multi-element ocean wave data, characterized in that, Includes the following steps: A multi-element wave state tensor space is constructed, wherein the multi-element wave state tensor space uses longitude, latitude, depth and time as coordinate axes and multiple wave elements as channels; The phase alignment process is performed on each element channel in the multi-element wave state tensor space by the tidal phase modulation function to generate a phase-aligned tensor. The phase-aligned tensor is input into a multi-factor assimilation model, which includes an encoder, a spectral decomposition layer, a coupled iterative layer, and a decoder. The encoder encodes the phase-aligned tensor to obtain latent features. The spectral decomposition layer decomposes the latent features channel by channel to obtain spectral decomposition features. The coupled iterative layer iterates the spectral decomposition features across channels to obtain updated latent features. The decoder decodes the updated latent features to obtain the assimilated wave element tensor. Based on the assimilated wave element tensor, a reanalysis field of wave elements for the target sea area is generated; The phase alignment processing of each element channel in the multi-element ocean wave state tensor space using the tidal phase modulation function is specifically as follows: Obtain the tidal harmonic constant of the target sea area; Construct a tidal phase field corresponding to each coordinate point in the multi-element wave state tensor space based on the tidal harmonic constant; For each element channel, calculate the phase difference between the instantaneous phase of that element channel and the tidal phase field; Construct a rotation matrix based on the phase difference; The rotation matrix is ​​used to rotate and project the feature channel in the depth and time dimensions to obtain phase-aligned channel data. All phase-aligned channel data are recombined along the channel dimension to generate the phase-aligned tensor.

2. The method for assimilation of multi-element ocean wave data according to claim 1, characterized in that, The construction of the multi-element wave state tensor space is specifically as follows: Acquire raw wave data from multiple observation sources for the target sea area. The wave data include significant wave height, mean wave direction, mean wave period, and sea surface wind speed. The observation sources include satellite remote sensing data, buoy observation data, and numerical model output data. The target sea area is divided into unstructured spatial grids to obtain a spatial grid, and the observation time range is discretized to establish a unified time series; The raw data of each wave element is mapped onto the spatial grid and the time series respectively to obtain the initial three-dimensional data field of each wave element; The initial three-dimensional data fields of all wave elements are stacked along the element dimensions to generate the multi-element wave state tensor space.

3. The method for assimilation of multi-element ocean wave data according to claim 2, characterized in that, The process of mapping the raw data of each wave element onto the spatial grid and the time series specifically involves: For the raw data of each wave element, select the grid points of the spatial grid; A neighborhood window is constructed centered on the grid points; Within the neighborhood window, the original data is subjected to a discrete step length weighted average to obtain the feature value at the grid point; Repeat the above weighted average operation for each time point in the time series to obtain the initial three-dimensional data field of the wave element.

4. The method for assimilation of multi-element ocean wave data according to claim 1, characterized in that, The step of inputting the phase-aligned tensor into the multi-factor assimilation model specifically involves: The phase alignment tensor is input into the encoder, and the encoder performs three-dimensional convolutional encoding on the phase alignment tensor to obtain the first latent feature. The first latent feature is input into the spectral decomposition layer, and the spectral decomposition layer performs channel-by-channel spectral decomposition on the first latent feature to obtain high-frequency components and low-frequency components. The high-frequency component and the low-frequency component are input into the coupling iteration layer. The coupling iteration layer performs cross-channel feature coupling iteration on the high-frequency component and the low-frequency component, and recombines the iterated high-frequency component and the iterated low-frequency component along the feature channel dimension to obtain the second potential feature. The second latent feature is input into the decoder, which performs deconvolution decoding on the second latent feature to obtain the assimilated wave element tensor.

5. The method for assimilation of multi-element ocean wave data according to claim 4, characterized in that, The encoder performs three-dimensional convolutional encoding on the phase-aligned tensor, specifically as follows: The phase alignment tensor is divided into blocks along four dimensions: longitude, latitude, depth, and time, to obtain multiple spatial-temporal blocks; Each spatial-temporal block is input into the first three-dimensional convolutional layer, which contains eight convolutional kernels. The size of each convolutional kernel is 3, 3, 3, and 2 along the longitude, latitude, depth, and time, respectively. Batch normalization and nonlinear activation are applied to the output of the first three-dimensional convolutional layer to obtain the first convolutional feature. The first convolutional feature is input into a second three-dimensional convolutional layer, which contains 16 convolutional kernels. The size of each convolutional kernel is 3, 3, 1, and 2 along the longitude, latitude, depth, and time, respectively. Max pooling is performed on the output of the second 3D convolutional layer, with the pooling window having a length, latitude, depth, and time of 2, 2, 1, and 2, respectively, to obtain the first latent feature.

6. The method for assimilation of multi-element ocean wave data according to claim 4, characterized in that, The spectral decomposition layer performs channel-by-channel spectral decomposition on the first latent feature to obtain high-frequency and low-frequency components, specifically: For each element channel of the first potential feature, perform a three-level wavelet packet decomposition along the time axis to obtain 8 frequency band sub-sequences; For each frequency band subsequence, the sum of squares of each sampling point in the frequency band subsequence is taken as the frequency band energy. Frequency band subsequences whose frequency band energy exceeds the average frequency band energy of the element channel are marked as high frequency subsequences, and the remaining frequency band subsequences are marked as low frequency subsequences. The frequency band subsequences marked as high-frequency subsequences are reconstructed by inverse wavelet packet to obtain the high-frequency tensor, and the frequency band subsequences marked as low-frequency subsequences are reconstructed by inverse wavelet packet to obtain the low-frequency tensor. The high-frequency tensor and the low-frequency tensor are concatenated according to the element channel dimension to obtain the high-frequency component and the low-frequency component.

7. The method for assimilation of multi-element ocean wave data according to claim 4, characterized in that, The coupling iteration layer performs cross-channel feature coupling iteration on the high-frequency components and the low-frequency components, specifically as follows: The high-frequency component and the low-frequency component are spliced ​​together along the element channel dimension to obtain the feature to be coupled. For any two feature channels in the features to be coupled, calculate the Pearson correlation coefficient between feature values ​​at the same spatial location and time point, and construct a cross-correlation matrix using the Pearson correlation coefficients between all channel pairs; A channel coupling graph is constructed based on the cross-correlation matrix, where the nodes of the channel coupling graph are feature channels and the edge weights are the Pearson correlation coefficients between the corresponding channels. Perform spectral clustering based on the normalized Laplacian matrix on the channel coupling graph to obtain multiple coupled channel groups; The channel features within each coupled channel group are weighted and averaged, and the weighted average result is used to update the features of all channels within that coupled channel group. The updated channel features are used again as the features to be coupled, and the cross-correlation matrix calculation, channel coupling graph construction, spectral clustering and weighted average update are performed iteratively until the preset number of iterations is reached. After the final iteration, the updated channel features are split into high-frequency components and low-frequency components according to the splicing order.

8. The method for assimilation of multi-element ocean wave data according to claim 1, characterized in that, The step of generating the wave element reanalysis field of the target sea area based on the assimilated wave element tensor specifically involves: The assimilated wave element tensor is split along the channel dimension to obtain the assimilated data volume of each wave element; Perform an inverse transform on the assimilated data volume of each wave element, corresponding to the channel-by-channel spectral decomposition performed by the spectral decomposition layer, to obtain the spatial distribution field of the wave element. Based on the coastline boundary and seabed topography data of the target sea area, the spatial distribution field is modified by boundary constraints to obtain a modified spatial distribution field. The boundary constraint modification includes setting the feature values ​​at the coastline boundary to zero flux boundary conditions and setting the feature values ​​at the seabed topography to topography-following boundary conditions. The corrected spatial distribution fields of all wave elements are combined and interpolated onto the regular latitude and longitude grid of the target sea area to generate the wave element reanalysis field.

9. A method for assimilating multi-element ocean wave data according to claim 8, characterized in that, The interpolation onto the regular latitude and longitude grid of the target sea area specifically involves: Obtain the discrete point coordinates and discrete point values ​​of the modified spatial distribution field; Using each grid point of the regular latitude and longitude grid as the center, discrete points within a horizontal search radius of 50 kilometers are used as the interpolation neighborhood. For discrete points within the interpolation neighborhood, the spatial distance from the discrete point to the grid point is used as the weight, with the smaller the spatial distance, the greater the weight. The weighted average of the discrete point values ​​is calculated to obtain the initial interpolation at the grid point. The initial interpolation is smoothed by Gaussian filtering to obtain smoothed grid point values. The kernel size of the Gaussian filter is 3×3 along both the longitude and latitude directions. The smoothed grid point values ​​are used as the values ​​of the wave element reanalysis field on the regular latitude and longitude grid.

Images