Co-track time-series sea surface wind speed retrieval method and system based on GNSS-R and scatterometer

Abstract

The present invention belongs to the technical field of ocean data processing and discloses a co-track time-series sea surface wind speed inversion method and system combined with GNSS-R and scatterometer. The method comprises the following steps: acquiring historical satellite data of ocean surface wind, preprocessing the historical satellite data to obtain preprocessed data, and constructing a training sample set in combination with reference wind speed and surge auxiliary data; constructing an initial wind speed inversion model by using a plurality of single-moment multi-frequency information fusion modules and a multi-moment interactive attention enhancement module, and training the initial wind speed inversion model by using the training sample set to obtain a co-track time-series sea surface wind speed inversion model; acquiring a delay-Doppler correlation power diagram, GNSS-R observation parameters, scatterometer parameters, and wind wave and surge significant wave height, and completing sea surface wind speed inversion by using the co-track time-series sea surface wind speed inversion model to obtain sea surface wind speed prediction data.

Description

Technical Field

[0001] The present invention belongs to the technical field of ocean data processing, and in particular relates to a method and system for inverting sea surface wind speed in a co-track time series manner using a combined GNSS-R and scatterometer. Background Art

[0002] The L-band signals used by GNSS satellites have strong penetration capabilities and are less affected by clouds, fog, rain, and snow, making them a free signal source for Earth observation. GNSS signals are affected by the Earth's surface during reflection, causing variations in signal strength, latency, and frequency. Ocean surface wind field is a key factor affecting sea surface roughness, with surface roughness varying under different wind speeds. These differences in roughness cause variations in reflectivity, which in turn alters the power and frequency of the reflected signal, ultimately capturing it in the receiver. GNSS-R wind speed inversion utilizes the relationship between wind field, roughness, reflectivity, and signal characteristics.

[0003] The current GNSS-R wind speed inversion can be divided into three categories: (1) Wind speed inversion method based on waveform matching: Waveform matching is one of the commonly used methods for GNSS-R wind speed inversion. The implementation process of this method is to first obtain a series of DDMs through field measurements or numerical simulations and record the corresponding wind speed values ​​to build a "DDM-wind speed" benchmark library for matching. Then, the DDM of the unknown wind speed is compared (matched) with each DDM in the benchmark library, and the wind speed corresponding to the benchmark DDM with the smallest difference is used as the inverted wind speed. However, this method requires the establishment of a detailed theoretical waveform library, which is computationally intensive, time-consuming, and the inversion process is cumbersome. (2) Wind speed inversion method based on geographic model function: The wind speed inversion method based on empirical function mainly inverts wind speed by establishing a mapping relationship between GNSS-R observations and wind speed. This method is also one of the commonly used methods for GNSS-R wind speed inversion. The empirical function used here is the Geography Model Function (GMF), which includes exponential functions, polynomial functions, and 2D lookup tables (LUTs). However, this method mainly relies on a single observation value in the DDM, and it is difficult to comprehensively consider the influence of the reflection geometric parameters. (3) Wind speed inversion method based on machine learning (deep learning): Compared with the empirical function model, the machine learning model has a stronger nonlinear mapping ability. This type of method uses DDM data and its derived features as the main input, and mines deep features through models such as random forest, BP neural network, and convolutional neural network to establish a mapping relationship between input and wind speed. Although this method has strong nonlinear fitting ability, it has the problems of insufficient single-source feature information and failure to integrate multi-source remote sensing data. In addition, it is mostly based on single-time modeling, ignoring the temporal characteristics of multi-time observations along the track, which limits the inversion accuracy and model generalization ability.

[0004] To solve the above problems, the present invention extracts multi-frequency information based on the DDM characteristics and fuses multi-source observation data, and combines the interactive attention mechanism for time series modeling, aiming to improve the accuracy and robustness of sea surface wind speed inversion. Summary of the Invention

[0005] The present invention aims to solve the deficiencies of the prior art and provides the following solutions:

[0006] The GNSS-R and scatterometer combined co-track time series sea surface wind speed inversion method includes the following steps:

[0007] Acquiring historical satellite data of ocean surface winds, preprocessing the historical satellite data to obtain preprocessed data, and constructing a training sample set in combination with reference wind speed and surge auxiliary data;

[0008] An initial wind speed inversion model is constructed using several single-moment multi-frequency information fusion modules and multi-moment interactive attention enhancement modules, and the initial wind speed inversion model is trained using the training sample set to obtain a co-track time series sea surface wind speed inversion model;

[0009] The delay-Doppler correlation power diagram, GNSS-R observation parameters, scatterometer parameters and wind wave and surge significant wave height are obtained, and the sea surface wind speed inversion is completed using the same-track time series sea surface wind speed inversion model to obtain sea surface wind speed prediction data.

[0010] Preferably, the method for preprocessing the historical satellite data includes:

[0011] The range gain correction is calculated by the antenna gain at the mirror reflection point, the distance from the receiver to the mirror reflection point, and the distance from the navigation satellite to the mirror reflection point, and the range gain correction is performed on the historical satellite data:

[0012] ,

[0013] Where RCG represents the range gain correction, G SP represents the antenna gain in the direction of the mirror emission point, Indicates the distance from the receiver to the mirror reflection point, R SP Indicates the distance from the navigation satellite to the mirror reflection point, 10 27 represents the scaling factor;

[0014] Performing incident angle correction on the corrected data, and filtering out data having an incident angle less than a preset threshold value in the corrected data;

[0015] Perform specular reflection point type correction on the filtered data. In the specular reflection point type, 0 represents open sea area, 0.5 represents coastal sea area within 25 kilometers of land, 1 represents land, and 2 represents sea ice. Filter out the data with the specular reflection point type less than 1 to obtain the type-corrected data;

[0016] A multi-look mean operation is performed on the wind field measurement radar data, backscatter coefficients of the wind field measurement radar data under C-band HH polarization and VV polarization are extracted, and the type-corrected data and the scatter coefficients are integrated to obtain the preprocessed data.

[0017] Preferably, the method for constructing the training sample set includes:

[0018] The U component and V component of the 10-meter wind field measured by ECMWF are selected, and the reference wind speed is calculated based on the U component and the V component:

[0019] ,

[0020] Among them, wspd represents the reference wind speed, u 10 represents the U component of the 10-meter wind field, v 10 represents the V component of the 10-meter wind field;

[0021] The significant wave heights of wind waves and swell waves are selected as the swell auxiliary data, and the reference wind speed and the swell auxiliary data are bilinearly interpolated to obtain the training sample set.

[0022] Preferably, the single-time multi-frequency information fusion module includes: a wavelet transform unit, a global enhancement unit, a local enhancement unit, a BP unit and a feature fusion unit;

[0023] The wavelet transform unit is composed of a low-pass filter and a high-pass filter; the wavelet transform unit is used to perform two-dimensional wavelet decomposition on the delay-Doppler correlation power map, divide the image into one low-frequency approximate sub-band and three high-frequency detail sub-bands, and use the low-frequency approximate sub-band as a low-frequency DDM image, and perform weighted fusion on the three high-frequency detail sub-bands to obtain a high-frequency DDM image;

[0024] The global enhancement unit is composed of tile embedding, position encoding, layer normalization, multi-head self-attention and multi-layer perceptron; the tile embedding is used to convert the low-frequency DDM image and the high-frequency DDM image into a feature vector sequence; the position encoding is used to introduce a delay Doppler coordinate with the same size as the feature vector sequence, and integrate the feature vector sequence with the delay Doppler coordinate to obtain an input vector; the layer normalization is used to normalize the input vector to obtain a normalized vector; the multi-head self-attention processes the normalized vector based on the self-attention mechanism to obtain a number of input features; the multi-layer perceptron is used to perform linear transformation and nonlinear activation on each of the input features to obtain a global feature sequence;

[0025] The local enhancement unit includes a CNN network and a bidirectional feature enhancement subunit; the CNN is used to preliminarily extract local texture features of the low-frequency DDM image and the high-frequency DDM image, and the bidirectional feature enhancement subunit is used to perform a global average pooling operation on the local texture features and weight them to obtain a local feature sequence;

[0026] The BP unit is used to extract features of GNSS-R observation parameters, scatterometer parameters and wind wave and surge significant wave height and perform weighted summation to obtain multi-source parameter fusion features;

[0027] The feature fusion unit is used to fuse the global feature sequence, the local feature sequence and the multi-source parameter fusion feature to obtain a single-moment fusion feature.

[0028] Preferably, the multi-moment interactive attention enhancement module includes: a first bidirectional long short-term memory network, a second bidirectional long short-term memory network and a fully connected layer;

[0029] The first bidirectional long short-term memory network is used to extract the time-dependent features of the single-moment fusion features, and combines the multi-head self-attention mechanism to enhance the internal correlation between different time steps;

[0030] The second bidirectional long short-term memory network is used to further extract the deep temporal features of the single-moment fusion features, and strengthen the interaction modeling between features at different time steps through a multi-head cross attention mechanism;

[0031] The fully connected layer is used to fuse the time-dependent features and the deep temporal features to obtain a sea surface wind speed prediction result.

[0032] Preferably, the method for training the initial wind speed inversion model includes:

[0033] Initializing parameters in the initial wind speed inversion model, and transferring the training sample set into the initial wind speed inversion model;

[0034] The training sample set obtains a predicted value through forward propagation and calculates a loss function;

[0035] Using the back propagation algorithm and the Adam optimizer to update the parameters in the initial wind speed inversion model;

[0036] When the change of the continuous iterative loss function is less than 0.0001 or the set number of iterations is reached, the training ends and the same-track multi-time series sea surface wind speed inversion model is obtained.

[0037] Preferably, the method for calculating the loss function includes: constructing the loss function using a mean square error loss function, a direction consistency loss, and an amplitude consistency loss:

[0038] ,

[0039] ,

[0040] ,

[0041] ,

[0042] Among them, L Total represents the loss function, L MSE represents the mean square error loss function, L Direction represents the direction consistency loss, L Trend represents the amplitude consistency loss, α, β, and γ represent weight coefficients, softplus represents the softplus function, w1, w2, and w3 represent learnable parameters;

[0043] The calculation method of the mean square error loss function includes:

[0044] ,

[0045] Among them, n represents the number of samples, y i represents the true value of the i-th sample, f i represents the inversion value of the i-th sample;

[0046] The calculation method of the directional consistency loss includes:

[0047] ,

[0048] ,

[0049] ,

[0050] Where Δy i Represents the difference between the true value of the sample, Δf iIndicates the difference of sample inversion values, y i+1 represents the true value of the i+1th sample, f i+1 represents the inversion value of the i+1th sample;

[0051] The calculation method of the amplitude consistency loss includes:

[0052] ,

[0053] ,

[0054] ,

[0055] Among them, Var(Δy) represents the variance of the true value of the sample, Var(Δf) represents the variance of the inverted value of the sample, represents the mean of the sample true value differences, Represents the mean of the sample inversion value differences.

[0056] Preferably, the update rules of the Adam optimizer include:

[0057] ,

[0058] Among them, t represents the number of iterations, g t represents the current gradient, Represents the gradient of the loss function L with respect to the model parameter 𝜔, m t represents the first-order moment estimate of the gradient, β1 represents the decay rate of the control first-order moment estimate, m t-1 represents the first-order moment estimate at the previous iteration, v t represents the second-order moment estimate of the gradient, β2 represents the decay rate of the control second-order moment estimate, It represents the result after correcting the estimated value of the first-order moment estimate. represents the decay rate β1 to the power of t, It represents the result after correcting the estimated value of the second-order moment estimate. represents the attenuation rate β1 to the power of t, Δω t represents the updated value of the model parameter 𝜔 at the tth iteration, α1 represents the learning rate, represents the smoothing term.

[0059] The present invention also provides a co-track time-series sea surface wind speed inversion system combined with GNSS-R and scatterometer, wherein the inversion system applies any of the above-mentioned inversion methods and comprises: a data set construction module, a model construction module and a prediction module;

[0060] The data set construction module is used to obtain historical satellite data of ocean surface wind, preprocess the historical satellite data to obtain preprocessed data, and construct a training sample set in combination with reference wind speed and surge auxiliary data;

[0061] The model construction module uses a plurality of single-time multi-frequency information fusion modules and a multi-time interactive attention enhancement module to construct an initial wind speed inversion model, and uses the training sample set to train the initial wind speed inversion model to obtain a co-track time series sea surface wind speed inversion model;

[0062] The prediction module is used to obtain the delay-Doppler correlation power diagram, GNSS-R observation parameters, scatterometer parameters and wind wave and surge effective wave height, and use the same track time series sea surface wind speed inversion model to complete the sea surface wind speed inversion to obtain sea surface wind speed prediction data.

[0063] Compared with the prior art, the present invention has the following beneficial effects:

[0064] Based on the structural characteristics of GNSS-R DDM images, this paper designs a single-moment multi-frequency information fusion module (MIFM) to extract complementary features between high- and low-frequency DDM images and multi-source physical parameters. Furthermore, considering the intra-orbital spatiotemporal continuity of spaceborne observation data, to enhance the model's ability to model the temporal evolution of wind speed, this paper further introduces self-attention and cross-attention mechanisms on the basis of the traditional BiLSTM architecture, constructing a multi-moment interactive attention enhancement module (MIAE). This strengthens the interactive modeling and dynamic association between temporal features. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] In order to more clearly illustrate the technical solution of the present invention, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0066] Figure 1 Schematic diagram of a method flow in an embodiment of the present invention;

[0067] Figure 2 Schematic diagram of the structure of the initial wind speed inversion model according to an embodiment of the present invention;

[0068] Figure 3 A single-time multi-frequency information fusion module according to an embodiment of the present invention;

[0069] Figure 4 This is a schematic diagram of the structure of a bidirectional feature enhancer unit according to an embodiment of the present invention;

[0070] Figure 5Schematic diagram of the BP unit structure according to an embodiment of the present invention;

[0071] Figure 6 This is a multi-moment interactive attention enhancement module according to an embodiment of the present invention. DETAILED DESCRIPTION

[0072] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0073] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0074] Example 1

[0075] In this embodiment, if Figure 1 As shown in FIG, the inversion method of the sea surface wind speed in the same track time series by combining GNSS-R and scatterometer includes the following steps:

[0076] S1. Obtain historical satellite data of ocean surface winds, preprocess the historical satellite data to obtain preprocessed data, and construct a training sample set based on the reference wind speed and surge auxiliary data.

[0077] In this embodiment, historical satellite data of ocean surface winds from Tianmu-1 (TM) are first obtained and preprocessed. The preprocessing method includes calculating a range gain correction based on the antenna gain at the mirror reflection point, the distance from the receiver to the mirror reflection point, and the distance from the navigation satellite to the mirror reflection point, and performing the range gain correction on the historical satellite data:

[0078] ,

[0079] Where RCG represents the range gain correction, G SP represents the antenna gain in the direction of the mirror emission point, Indicates the distance from the receiver to the mirror reflection point, R SP Indicates the distance from the navigation satellite to the mirror reflection point, 10 27Represents the scaling factor, used to adjust the RCG to a range of 1-100. To eliminate data with low signal gain, this paper selects data with an RCG greater than 10 for subsequent wind speed inversion to reduce the interference of low-quality data on the inversion results. The corrected data is corrected for the angle of incidence, and data with an angle of incidence less than a preset threshold is selected. In this embodiment, since antenna gain decreases with increasing angle of incidence, and the quality of the reflected signal decreases accordingly, data with an angle of incidence less than 60° is selected during data processing to ensure signal quality. The filtered data were corrected for the mirror reflection point type. In the mirror reflection point type, 0 represents open sea area, 0.5 represents coastal sea area within 25 kilometers from land, 1 represents land, and 2 represents sea ice. The data with the mirror reflection point type less than 1 were filtered out to obtain the type-corrected data. The multi-look mean operation was performed on the wind field measurement radar data, and the backscattering coefficient of the wind field measurement radar data under C-band HH polarization and VV polarization was extracted. The type-corrected data and the scattering coefficient were integrated to obtain the preprocessed data.

[0080] To ensure consistency and matching accuracy in multi-source data fusion, this embodiment uses a 15-minute time window for synchronization processing. The spatiotemporal matching conditions for TM data and FY data from the FY-3E wind radar are: a latitude and longitude difference of less than 0.1 degrees and a time difference of less than 15 minutes. After matching, the backscatter coefficient provided by FY is interpolated to the corresponding TM data point. Given that a single TM data point may correspond to multiple FY data points, and that closer observations have a higher weight in the interpolation result, this embodiment uses the inverse distance weighted (IDW) interpolation method to weight the FY data. Ultimately, an estimate of the wind radar backscatter coefficient at each TM data specular reflection point is obtained.

[0081] The method of constructing the training sample set includes: the reference wind speed is the east-west wind speed component U at a height of 10 meters provided by ECMWF 10 (10m u-component of wind) and the north-south wind speed component V at a height of 10 meters 10 (10m v-component of wind) is calculated, and the wind speed formula is as follows:

[0082] ,

[0083] Among them, wspd represents the reference wind speed, u 10 represents the U component of the 10-meter wind field, v 10 Represents the V component of the 10-meter wind field.

[0084] The significant wave height of wind waves and swells is selected as auxiliary swell data, and the reference wind speed and swell auxiliary data are bilinearly interpolated to obtain a training sample set. Considering that the reference wind speed and SWH are both global grid data, this embodiment uses a bilinear interpolation method to map their values ​​to the location of the TM data mirror reflection point. The specific calculation process is as follows: To obtain the estimated value of the unknown function f at P(x, y), it is known that Q 11 (x1,y1),Q 12 (x1,y2),Q 21 (x2,y1),Q 22 The coordinates of (x2, y2) and the corresponding pixel values. First, interpolate in the x direction to get R1 and R2:

[0085] ,

[0086] ,

[0087] Among them, R1 represents Q 11 and Q 21 The point obtained by horizontal interpolation, R2 represents Q 12 and Q 22 The points obtained by horizontal interpolation, f(R1) represents the function value at point R1, f(R2) represents the function value at point R2, f(Q 11 ) means at point Q 11 The function value at f(Q 12 ) means at point Q 12 The function value at f(Q 21 ) means at point Q 21 The function value at f(Q 22 ) means at point Q 22 The function value at , and then interpolate in the y direction based on R1 and R2 to obtain the function value f(P) of the target point P to be interpolated:

[0088] ,

[0089] In this embodiment, the target point P is the TM longitude and latitude coordinate data point of the mirror reflection point of the TM; the reference wind speed data and the SWH data are bilinearly interpolated to the TM data point to complete the data set creation.

[0090] S2. Use several single-time multi-frequency information fusion modules and multi-time interactive attention enhancement modules to build the initial wind speed inversion model, such as Figure 2 As shown in the figure, the initial wind speed inversion model is trained using the training sample set to obtain the same-track time series sea surface wind speed inversion model.

[0091] Single-time multi-frequency information fusion module Figure 3As shown, it includes: wavelet transform unit, global enhancement unit, local enhancement unit, BP unit and feature fusion unit.

[0092] The wavelet transform unit consists of a low-pass filter and a high-pass filter. The wavelet transform unit is used to perform two-dimensional wavelet decomposition on the delay-Doppler correlation power map, dividing the image into one low-frequency approximate subband and three high-frequency detail subbands. The low-frequency approximate subband is used as the low-frequency DDM image, and the three high-frequency detail subbands are weightedly fused to obtain the high-frequency DDM image.

[0093] In this embodiment, in order to extract different scale features in the DDM image, this embodiment performs wavelet transform on the DDM image and decomposes it into high-frequency and low-frequency components using a two-dimensional discrete wavelet transform. Specifically, this embodiment uses Haar wavelet to perform a two-dimensional wavelet decomposition on the input DDM image, dividing the image into one low-frequency approximate subband and three high-frequency detail subbands, which are as follows: Low-frequency approximate subband:

[0094] ,

[0095] Horizontal high-frequency detail subband:

[0096] ,

[0097] Vertical high-frequency detail subband:

[0098] ,

[0099] Diagonal high-frequency detail subband:

[0100] ,

[0101] Where X[k,l] is the pixel value of the original image in the 𝑘th row and the lth column, g[·] is the coefficient of the Haar wavelet low-pass filter, h[·] is the coefficient of the Haar wavelet high-pass filter, and [m,n] is the position index of the output sub-image; 2m and 2n indicate that the downsampling factor is 2. The filter coefficients of the Haar wavelet are defined as follows: Low-pass filter:

[0102] ,

[0103] High-pass filter:

[0104] ,

[0105] The low-frequency approximate subband cA obtained by wavelet transform is the low-frequency DDM image; the detail information in the three high-frequency subbands is formed into a high-frequency DDM image through weighted fusion. The fusion calculation formula is as follows:

[0106] ,

[0107] Taking the DDM image with an input size of 20×9 as an example, after two-dimensional wavelet transform processing, the low-frequency DDM and high-frequency DDM with a size of 10×4 are obtained respectively.

[0108] The global enhancement unit consists of tile embedding, position encoding, layer normalization, multi-head self-attention and multi-layer perceptron; tile embedding is used to convert low-frequency DDM images and high-frequency DDM images into feature vector sequences; position encoding is used to introduce delay Doppler coordinates that are consistent with the size of the feature vector sequence, and integrate the feature vector sequence with the delay Doppler coordinates to obtain the input vector; layer normalization is used to standardize the input vector to obtain a normalized vector; multi-head self-attention processes the normalized vector based on the self-attention mechanism to obtain several input features; multi-layer perceptron is used to perform linear transformation and nonlinear activation on each input feature to obtain a global feature sequence.

[0109] In this embodiment, the global enhancement unit establishes information connections between various positions of the image through the self-attention mechanism, so as to better extract the global features of high- and low-frequency DDM images. The global enhancement module designed in this embodiment consists of patch embedding, position encoding, layer normalization, multi-head self-attention and multi-layer perceptron. The details of each part will be introduced in detail below. The purpose of patch embedding is to convert image features into sequence input. In this module, the input image is first divided into blocks. In order to accurately capture the information of each pixel, the patch size is set to 1×1, and pixel-by-pixel encoding is performed to ensure that the local features of all pixels are fully preserved. For high-frequency or low-frequency DDM images after wavelet transformation, the size is 10×4, and the corresponding total number of patches is N=H×W=10×4=40. Each patch is recorded as patch i , where i=1,2,…,N. Flatten each patch to form a vector x i :

[0110] ,

[0111] Combine all flattened vectors into a feature vector sequence X:

[0112] ,

[0113] The main function of position encoding is to provide the model with the position information of the elements in the sequence. To enhance the model's ability to perceive global information, the delay-Doppler coordinate P with the same size as X is introduced for position marking:

[0114] ,

[0115] Ultimately, the encoder’s input vector Z is the sum of the feature vector X and the position marker P:

[0116] .

[0117] Layer normalization is used to standardize feature values ​​to a relatively stable range to accelerate the convergence process of model training.

[0118] The multi-head self-attention mechanism aims to fully model the similarity between input vectors, thereby capturing global dependencies and improving the network's ability to perceive image context information. The specific working principle of the multi-head self-attention is as follows: First, use the linear transformation matrix W Q 、W K and W V Multiplying by the input vector Z yields Q (Query), K (Key), and V (Value):

[0119] ,

[0120] ,

[0121] ,

[0122] Among them, Q represents the query vector, which is used to determine the information that should be focused on in the attention mechanism; K represents the key vector, which is used to calculate the correlation between Q and other position features in Z and provide additional context information. V represents the value vector, which contains the position information of each element in Z. Q 、W K and W V Randomly generated at the beginning of network training and updated during training:

[0123] ,

[0124] Among them, learning_rate represents the learning rate, Represents the derivative of the loss function with respect to the weight matrix; in the global enhancement unit, Q, K, and V are the key elements to implement the self-attention mechanism. The output of the self-attention mechanism is:

[0125] ,

[0126] Among them, Attention(Q,K,V) represents the attention weight, d k Indicates the dimension of the input vector. The essence of the multi-head self-attention mechanism is the linear transformation after the concatenation of multiple attention calculation results. This mechanism enables the model to use different feature information obtained at different positions, thereby increasing feature diversity. The calculation method of multi-head self-attention is as follows:

[0127] ,

[0128] ,

[0129] Where r represents the total number of heads. In this embodiment, the number of heads is 4. j Represents the features of the j-th head, W represents the corresponding weight matrix when obtaining the output result, and Concat represents the concatenation operation of the vectors.

[0130] To further enhance the module's nonlinear expressiveness, this example introduces a multilayer perceptron (MLP) after the attention output. This perceptron consists of two fully connected layers, with 128 and 64 neurons, respectively. It performs linear transformation and nonlinear activation on each input feature, thereby enhancing the model's expressiveness. Finally, after processing the high- and low-frequency DDM images through the global enhancement module, each outputs a global feature sequence with a feature length of 64, which serves as input for subsequent modules.

[0131] The local enhancement unit includes a CNN network and a bidirectional feature enhancement subunit; the CNN is used to preliminarily extract the local texture features of the low-frequency DDM image and the high-frequency DDM image, and the bidirectional feature enhancement subunit is used to perform a global average pooling operation on the local texture features and weight them to obtain a local feature sequence.

[0132] In this embodiment, the local enhancement unit extracts basic local texture features through CNN, and introduces a bidirectional feature enhancement module based on this. Through global average pooling in the row and column directions, it constructs attention weights in the horizontal and vertical directions, thereby enhancing the directional feature expression ability of the local area. The bidirectional feature enhancement module designed in this embodiment consists of global average pooling in the row and column directions, an attention weighting mechanism, and a feature flattening structure.

[0133] CNN is used to initially extract local texture features from high-frequency or low-frequency DDM images. The convolutional layer uses a 3×3 filter for feature extraction. This filter can be considered a local structure-aware feature extractor. It generates a feature map through convolution with the input image, which serves as the input for the bidirectional feature enhancement module. For example, for a high-frequency or low-frequency DDM image (4×10) as an input, the convolutional layer outputs a feature map of the same size (4×10×1). The filter output can be expressed as follows:

[0134] ,

[0135] Wherein, z represents the output result of the filter, T represents the input high-frequency or low-frequency DDM image, M and b represent the weight and bias of the convolution kernel, * represents the convolution operation, and f(·) represents the activation function, which is ReLU in this embodiment.

[0136] Bidirectional characteristic enhancer units such as Figure 4 As shown in the figure, it aims to improve the expression ability of the direction-sensitive significant area in the convolution feature map. Its input is the feature map after the above convolution processing, denoted as , where H=4 represents the height (number of rows) of the feature map F, and L=10 represents the width (number of columns) of the feature map F. First, global average pooling operations are applied in the vertical (row) and horizontal (column) directions to obtain two one-dimensional attention features, which are expressed as:

[0137] ,

[0138] ,

[0139] Multiply the two one-dimensional attention features through the outer product to construct a two-dimensional attention weight map:

[0140] ,

[0141] Where σ(·) represents the Sigmoid activation function, which is used to normalize the attention weight. Then, the attention weight is applied to the input feature map to obtain the weighted feature map:

[0142] ,

[0143] Finally, the weighted feature map is flattened into a local feature sequence with a feature length of 40 as the output of the local enhancement module.

[0144] The BP unit is used to extract features of GNSS-R observation parameters, scatterometer parameters and wind wave and surge significant wave height and perform weighted summation to obtain multi-source parameter fusion features.

[0145] In this embodiment, the BP unit is as follows Figure 5As shown in the figure, it is used for feature extraction of multi-source parameter information, and contains three fully connected layers. The input of this unit includes 9 categories of GNSS-R observation parameters, microwave scatterometer parameters and surge auxiliary information, namely the latitude and longitude of the mirror reflection point, normalized bistatic radar cross section (NBRCS), DDM leading edge slope (LES), signal-to-noise ratio, incident angle, range correction gain (RCG), HH polarization backscatter coefficient, VV polarization backscatter coefficient and significant wave height (SWH) of wind waves and surges. The three fully connected layers are composed of 32, 64 and 128 neurons respectively. Each neuron performs weighted summation with multiple output features of the previous layer and is processed by the activation function to generate an output feature. Therefore, the feature output by the i-th neuron of the first hidden layer can be expressed as:

[0146] ,

[0147] Among them, x j represents the j-th input feature of the neuron, Represents the weight value of the line connecting the neuron and the j-th input feature, represents the bias value of the neuron, and f(·) represents the activation function. In this embodiment, the activation function is ReLU.

[0148] The ReLU function is one of the most commonly used nonlinear functions in neural networks due to its advantages such as simple calculation, fast convergence speed, and effective mitigation of the vanishing gradient problem. This function can truncate all negative values ​​of the input data to 0 and output the positive values ​​intact. Its mathematical expression is:

[0149] ,

[0150] Among them, x represents the input data, and f(x) represents the output result of the activation function.

[0151] The feature fusion unit is used to fuse the global feature sequence, the local feature sequence and the multi-source parameter fusion features to obtain the single-moment fusion feature.

[0152] In this embodiment, the feature fusion unit concatenates the high-frequency and low-frequency DDM features extracted by the global enhancement module and the local enhancement module respectively to form a complete DDM spatial feature; then, it is concatenated again with the multi-source parameter fusion feature extracted by the BP module to construct a single-moment fusion feature.

[0153] Multi-moment interactive attention enhancement module such as Figure 6 As shown, it includes: a first bidirectional long short-term memory network, a second bidirectional long short-term memory network and a fully connected layer.

[0154] The first bidirectional long short-term memory network is used to extract time-dependent features from several single-moment fused features and, combined with a multi-head self-attention mechanism, enhance the internal correlations between different time steps. The second bidirectional long short-term memory network is used to further extract deep temporal features from the single-moment fused features and, using a multi-head cross-attention mechanism, strengthen the interaction modeling between features at different time steps. The fully connected layer is used to fuse the time-dependent features with the deep temporal features to produce the sea surface wind speed forecast.

[0155] In this embodiment, the core concept of BiLSTM is to introduce a backpropagation path based on the traditional LSTM, enabling the model to simultaneously utilize past and future sequence information for modeling, thereby enhancing its ability to model contextual dependencies. This bidirectional structure is particularly suitable for time series tasks that require comprehensive consideration of both historical and future contextual information. Similar to the traditional LSTM, BiLSTM uses the cell state as the primary carrier of information flow and controls the information update process through three gating mechanisms: input gate, forget gate, and output gate. These gated units generate weights between 0 and 1 using a sigmoid function to regulate the degree of acceptance of new information, the degree of retention of historical information, and the influence of the current state on the output. The synergistic effect of bidirectional information flow and gating structure enables BiLSTM to effectively capture sequence dependencies over long time spans, improving the model's performance in time series modeling tasks.

[0156] (a) Forget Gate: The forget gate is used to control the historical information components in the cell state that should be discarded. It generates a weight value between 0 and 1 through a sigmoid function, which represents the retention ratio of each state unit. Its mathematical expression is:

[0157] ,

[0158] Among them, f t represents the forget gate, σ represents the sigmoid function, W f Represents the weight matrix of the forget gate, h t-1 represents the hidden state of the previous time step, x t represents the input of the current time step, t represents the time step, b f Represents the bias vector of the forget gate.

[0159] (b) Input Gate: The input gate consists of two parts: a sigmoid layer, which determines the information component in the state that needs to be updated; and a tanh layer, which generates the candidate state vector at the current moment. The outputs of these two layers are element-wise multiplied to form the candidate information that is ultimately written into the cell state. The mathematical expression of the sigmoid layer is:

[0160] ,

[0161] The mathematical expression of the tanh layer is:

[0162] ,

[0163] Among them, i t represents the first sigmoid layer, W f represents the weight matrix of the first sigmoid layer, b f represents the bias vector of the first sigmoid layer, represents the first tanh layer, W C represents the weight matrix of the first tanh layer, b C Represents the bias vector of the first tanh layer.

[0164] (c) Cell State: The cell state is the core structure of LSTM, used to carry and store long-term dependency information across time steps. The cell state is updated by adding the output of the forget gate and the output of the input gate. Its mathematical expression is:

[0165] ,

[0166] Among them, C t Indicates the updated unit state, C t-1 Represents the cell state at the previous time step.

[0167] (d) Output Gate: The output gate controls the hidden layer output at the current time step. It uses a sigmoid layer to determine the state information to be included in the output, then forms a candidate output through a tanh layer, and finally combines these two parts to form the final output. The mathematical expression of the sigmoid layer of the output gate is:

[0168] ,

[0169] The mathematical expression of the final hidden layer state is:

[0170] ,

[0171] Among them, t represents the second sigmoid layer, Wo represents the weight matrix of the second sigmoid layer, b o Represents the bias vector of the second sigmoid layer, h t represents the final hidden layer.

[0172] The Multi-Moment Interactive Attention Enhancement (MIAE) module uses a two-layer BiLSTM network to model the temporal dependencies of multi-moment feature sequences and integrates self-attention and cross-attention mechanisms to enhance dynamic correlations and information interactions between different time steps. The module first receives the fused features of multiple time steps as a sequence input and feeds it into the first BiLSTM layer. It extracts a bidirectional latent state vector (128 dimensions) for each time step and constructs preliminary time-series-aware features. Subsequently, a self-attention mechanism is introduced to dynamically adjust the response weights of the features at each time step, highlighting the dominant role of key moments in wind speed fluctuations. The enhanced feature sequence is passed to the second BiLSTM layer to extract a deeper temporal feature representation (64 dimensions). Next, the module employs a cross-attention mechanism, using the output of the second BiLSTM layer as the query vector (Query) and the time-series-aware features of the first layer as the key (Key) and value (Value), achieving cross-modeling and feature fusion between different layers. Finally, the fused features are fed into a fully connected layer to output the sea surface wind speed retrieval results for consecutive time steps.

[0173] After the neural network model is constructed, it needs to be trained with a large number of historical samples to fit the nonlinear mapping relationship between input and output. The entire training process can be divided into two stages: forward propagation and backpropagation. Forward propagation involves the input data being weighted and processed by activation functions at each layer, and then passed layer by layer to the output layer to calculate the predicted value. Backpropagation, based on gradient descent and optimization algorithms, propagates the prediction error backward along the network structure and is used to update the parameter weights to minimize the overall loss.

[0174] During training, a loss function measures the difference between the model's inversion results and the true observations. The mean squared error (MSE) loss function is a commonly used evaluation metric in wind speed inversion models. However, in this embodiment, the model output is a sequence of wind speeds over multiple time steps. Using only the MSE metric makes it difficult to fully characterize the dynamic evolution of wind speed over time.

[0175] In this embodiment, the method for training the initial wind speed inversion model includes: initializing the parameters in the initial wind speed inversion model and passing the training sample set into the initial wind speed inversion model; obtaining a predicted value of the training sample set through forward propagation and calculating the loss function; using the backpropagation algorithm and the Adam optimizer to update the parameters in the initial wind speed inversion model; when the change in the continuous iterative loss function is less than 0.0001 or reaches the set number of iterations, the training ends and a same-track multi-time series sea surface wind speed inversion model is obtained.

[0176] The method of calculating the loss function includes: constructing the loss function using the mean square error loss function, direction consistency loss and amplitude consistency loss:

[0177] ,

[0178] ,

[0179] ,

[0180] ,

[0181] Among them, L Total represents the loss function, L MSE represents the mean square error loss function, L Direction represents the direction consistency loss, L Trend represents the amplitude consistency loss, α, β, and γ represent weight coefficients, softplus represents the softplus function, w1, w2, and w3 represent learnable parameters; the calculation method of the mean square error loss function includes:

[0182] ,

[0183] Among them, n represents the number of samples, y i represents the true value of the i-th sample, f i Represents the inversion value of the i-th sample; the calculation method of directional consistency loss includes firstly performing the difference between the predicted wind speed sequence and the true wind speed sequence, calculating the gradient change between adjacent time steps, forming the predicted gradient vector and the true gradient vector, and then using the cosine similarity to quantify the degree of directional similarity between the two sets of gradient vectors. The calculation formula is as follows:

[0184] ,

[0185] ,

[0186] ,

[0187] Where Δy iRepresents the difference between the true value of the sample, Δf i Indicates the difference of sample inversion values, y i+1 represents the true value of the i+1th sample, f i+1 represents the inversion value of the i+1th sample; the calculation method of the amplitude consistency loss includes calculating the variance of the predicted gradient and the true gradient, and taking the absolute difference between the two as the amplitude consistency loss, thereby prompting the model to be consistent with the actual observation in the local fluctuation amplitude. The calculation formula is as follows:

[0188] ,

[0189] ,

[0190] ,

[0191] Among them, Var(Δy) represents the variance of the true value of the sample, Var(Δf) represents the variance of the inverted value of the sample, represents the mean of the sample true value differences, Represents the mean of the sample inversion value differences.

[0192] This embodiment uses the Adaptive Moment Estimation (Adam) optimization algorithm, which is used to iteratively update the parameters in the neural network during the error backpropagation process to accelerate convergence and minimize the loss function. The update rules of the Adam optimizer are as follows:

[0193] ,

[0194] Among them, t represents the number of iterations, g t represents the current gradient, Represents the gradient of the loss function L with respect to the model parameter 𝜔, m t represents the first-order moment estimate of the gradient, β1 represents the decay rate of the control first-order moment estimate, m t-1 represents the first-order moment estimate at the previous iteration, v t represents the second-order moment estimate of the gradient, β2 represents the decay rate of the control second-order moment estimate, It represents the result after correcting the estimated value of the first-order moment estimate. represents the decay rate β1 to the power of t, It represents the result after correcting the estimated value of the second-order moment estimate. represents the attenuation rate β1 to the power of t, Δω t represents the updated value of the model parameter 𝜔 at the tth iteration, α1 represents the learning rate, represents the smoothing term.

[0195] S3. Obtain the delay-Doppler correlation power diagram, GNSS-R observation parameters, scatterometer parameters, and wind surge significant wave height, and use the same-track time-series sea surface wind speed inversion model to complete the sea surface wind speed inversion and obtain the sea surface wind speed forecast data.

[0196] In this embodiment, after completing the sea surface wind speed forecast, an appropriate regression evaluation metric is introduced to quantify the accuracy and fitting effect of the model inversion results. This embodiment uses the root mean square error (RMSE) as the primary evaluation metric. RMSE can be used to measure the overall deviation between the model prediction value and the reference value. Its mathematical expression is as follows:

[0197] ,

[0198] Where y i represents the true value of the i-th sample, f i represents the inversion value of the i-th sample; n represents the number of samples. A smaller RMSE value indicates a smaller difference between the predicted data and the reference data, and a better fitting effect; conversely, a larger difference between the predicted data and the reference data indicates a worse fitting effect.

[0199] Experiments were conducted based on data from September 2023, with the training and test sets divided in a ratio of 7:3. The time series length of a single input was set to N = 10. The experimental results show that the root mean square error (RMSE) of the proposed method is 1.13 m / s, demonstrating excellent performance in both prediction accuracy and time series modeling capabilities. The experimental results are shown in Table 1:

[0200] Table 1

[0201] .

[0202] Example 2

[0203] In this embodiment, a co-track time-series sea surface wind speed inversion system combining GNSS-R and scatterometer includes: a data set construction module, a model construction module and a prediction module.

[0204] The dataset construction module is used to obtain historical satellite data of ocean surface wind, preprocess the historical satellite data to obtain preprocessed data, and construct a training sample set in combination with reference wind speed and surge auxiliary data.

[0205] The model construction module uses several single-moment multi-frequency information fusion modules and multi-moment interactive attention enhancement modules to construct an initial wind speed inversion model, and uses the training sample set to train the initial wind speed inversion model to obtain the same-track time series sea surface wind speed inversion model.

[0206] The prediction module is used to obtain the delay-Doppler correlation power diagram, GNSS-R observation parameters, scatterometer parameters and wind wave and surge effective wave height, and use the same-track time-series sea surface wind speed inversion model to complete the sea surface wind speed inversion and obtain the sea surface wind speed prediction data.

[0207] The embodiments described above are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by persons skilled in the art should fall within the scope of protection defined by the claims of the present invention.

Claims

1. The co-track time series sea surface wind speed inversion method using GNSS-R and scatterometer is characterized by: The following steps are involved: Acquiring historical satellite data of ocean surface winds, preprocessing the historical satellite data to obtain preprocessed data, and constructing a training sample set in combination with reference wind speed and surge auxiliary data; An initial wind speed inversion model is constructed using several single-moment multi-frequency information fusion modules and multi-moment interactive attention enhancement modules, and the initial wind speed inversion model is trained using the training sample set to obtain a co-track time series sea surface wind speed inversion model; Obtaining the delay-Doppler correlation power diagram, GNSS-R observation parameters, scatterometer parameters, and wind wave and surge significant wave height, and using the same-track time series sea surface wind speed inversion model to complete the sea surface wind speed inversion and obtain sea surface wind speed prediction data; The single-time multi-frequency information fusion module includes: a wavelet transform unit, a global enhancement unit, a local enhancement unit, a BP unit and a feature fusion unit; The wavelet transform unit is composed of a low-pass filter and a high-pass filter; the wavelet transform unit is used to perform two-dimensional wavelet decomposition on the delay-Doppler correlation power map, divide the image into one low-frequency approximate sub-band and three high-frequency detail sub-bands, and use the low-frequency approximate sub-band as a low-frequency DDM image, and perform weighted fusion on the three high-frequency detail sub-bands to obtain a high-frequency DDM image; The global enhancement unit is composed of tile embedding, position encoding, layer normalization, multi-head self-attention and multi-layer perceptron; the tile embedding is used to convert the low-frequency DDM image and the high-frequency DDM image into a feature vector sequence; the position encoding is used to introduce a delay Doppler coordinate with the same size as the feature vector sequence, and integrate the feature vector sequence with the delay Doppler coordinate to obtain an input vector; the layer normalization is used to normalize the input vector to obtain a normalized vector; the multi-head self-attention processes the normalized vector based on the self-attention mechanism to obtain a number of input features; the multi-layer perceptron is used to perform linear transformation and nonlinear activation on each of the input features to obtain a global feature sequence; The local enhancement unit includes a CNN network and a bidirectional feature enhancement subunit; the CNN is used to preliminarily extract local texture features of the low-frequency DDM image and the high-frequency DDM image, and the bidirectional feature enhancement subunit is used to perform a global average pooling operation on the local texture features and weight them to obtain a local feature sequence; The BP unit is used to extract features of GNSS-R observation parameters, scatterometer parameters and wind wave and surge significant wave height and perform weighted summation to obtain multi-source parameter fusion features; The feature fusion unit is used to fuse the global feature sequence, the local feature sequence and the multi-source parameter fusion feature to obtain a single-moment fusion feature; The multi-moment interactive attention enhancement module includes: a first bidirectional long short-term memory network, a second bidirectional long short-term memory network and a fully connected layer; The first bidirectional long short-term memory network is used to extract the time-dependent features of the single-moment fusion features, and combines the multi-head self-attention mechanism to enhance the internal correlation between different time steps; The second bidirectional long short-term memory network is used to further extract the deep temporal features of the single-moment fusion features, and strengthen the interaction modeling between features at different time steps through a multi-head cross attention mechanism; The fully connected layer is used to fuse the time-dependent features and the deep temporal features to obtain a sea surface wind speed prediction result.

2. The method for inverting sea surface wind speed in a co-track time series by combining GNSS-R and scatterometer according to claim 1 is characterized in that: The method for preprocessing the historical satellite data includes: The range gain correction is calculated by the antenna gain at the mirror reflection point, the distance from the receiver to the mirror reflection point, and the distance from the navigation satellite to the mirror reflection point, and the range gain correction is performed on the historical satellite data: ， Where RCG represents the range gain correction, G SP represents the antenna gain in the direction of the mirror emission point, Indicates the distance from the receiver to the mirror reflection point, R SP Indicates the distance from the navigation satellite to the mirror reflection point, 10 27 represents the scaling factor; Performing incident angle correction on the corrected data, and filtering out data having an incident angle less than a preset threshold value in the corrected data; Perform specular reflection point type correction on the filtered data. In the specular reflection point type, 0 represents open sea area, 0.5 represents coastal sea area within 25 kilometers of land, 1 represents land, and 2 represents sea ice. Filter out the data with the specular reflection point type less than 1 to obtain the type-corrected data; A multi-look mean operation is performed on the wind field measurement radar data, backscatter coefficients of the wind field measurement radar data under C-band HH polarization and VV polarization are extracted, and the type-corrected data and the scatter coefficients are integrated to obtain the preprocessed data.

3. The method for inverting sea surface wind speed in a co-track time series by combining GNSS-R and scatterometer according to claim 1 is characterized in that: The method for constructing the training sample set includes: The U component and V component of the 10-meter wind field measured by ECMWF are selected, and the reference wind speed is calculated based on the U component and the V component: ， Among them, wspd represents the reference wind speed, u 10 represents the U component of the 10-meter wind field, v 10 represents the V component of the 10-meter wind field; The significant wave heights of wind waves and swell waves are selected as the swell auxiliary data, and the reference wind speed and the swell auxiliary data are bilinearly interpolated to obtain the training sample set.

4. The method for inverting sea surface wind speed in a co-track time series by combining GNSS-R and scatterometer according to claim 1 is characterized in that: The method for training the initial wind speed inversion model includes: Initializing parameters in the initial wind speed inversion model, and transferring the training sample set into the initial wind speed inversion model; The training sample set obtains a predicted value through forward propagation and calculates a loss function; Using the back propagation algorithm and the Adam optimizer to update the parameters in the initial wind speed inversion model; When the change of the continuous iterative loss function is less than 0.0001 or the set number of iterations is reached, the training ends and the same-track multi-time series sea surface wind speed inversion model is obtained.

5. The method for inverting sea surface wind speed in the same track time series by combining GNSS-R and scatterometer according to claim 4 is characterized in that: The method for calculating the loss function includes: constructing the loss function using a mean square error loss function, a direction consistency loss, and an amplitude consistency loss: ， ， ， ， Among them, L Total represents the loss function, L MSE represents the mean square error loss function, L Direction represents the direction consistency loss, L Trend represents the amplitude consistency loss, α, β, and γ represent weight coefficients, softplus represents the softplus function, w1, w2, and w3 represent learnable parameters; The calculation method of the mean square error loss function includes: ， Among them, n represents the number of samples, y i represents the true value of the i-th sample, f i represents the inversion value of the i-th sample; The calculation method of the directional consistency loss includes: ， ， ， Where Δy i Represents the difference between the true value of the sample, Δf i Indicates the difference of sample inversion values, y i+1 represents the true value of the i+1th sample, f i+1 represents the inversion value of the i+1th sample; The calculation method of the amplitude consistency loss includes: ， ， ， Among them, Var(Δy) represents the variance of the true value of the sample, Var(Δf) represents the variance of the inverted value of the sample, represents the mean of the sample true value differences, Represents the mean of the sample inversion value differences.

6. The method for inverting sea surface wind speed in a co-track time series by combining GNSS-R and scatterometer according to claim 4 is characterized in that: The update rules of the Adam optimizer include: ， Among them, t represents the number of iterations, g t represents the current gradient, Represents the gradient of the loss function L with respect to the model parameter 𝜔, m t represents the first-order moment estimate of the gradient, β1 represents the decay rate of the control first-order moment estimate, m t-1 represents the first-order moment estimate at the previous iteration, v t represents the second-order moment estimate of the gradient, β2 represents the decay rate of the control second-order moment estimate, It represents the result after correcting the estimated value of the first-order moment estimate. represents the decay rate β1 to the power of t, It represents the result after correcting the estimated value of the second-order moment estimate. represents the decay rate β1 to the power of t, Δω t represents the updated value of the model parameter 𝜔 at the tth iteration, α1 represents the learning rate, represents the smoothing term.

7. A co-track time-series sea surface wind speed inversion system combining GNSS-R and scatterometer, wherein the inversion system applies the inversion method according to any one of claims 1 to 6, characterized in that: include: Dataset building module, model building module and prediction module; The data set construction module is used to obtain historical satellite data of ocean surface wind, preprocess the historical satellite data to obtain preprocessed data, and construct a training sample set in combination with reference wind speed and surge auxiliary data; The model construction module uses a plurality of single-time multi-frequency information fusion modules and a multi-time interactive attention enhancement module to construct an initial wind speed inversion model, and uses the training sample set to train the initial wind speed inversion model to obtain a co-track time series sea surface wind speed inversion model; The prediction module is used to obtain the delay-Doppler correlation power diagram, GNSS-R observation parameters, scatterometer parameters and wind wave and surge effective wave height, and use the same track time series sea surface wind speed inversion model to complete the sea surface wind speed inversion to obtain sea surface wind speed prediction data.

Images