Method and System for Inverting Sea Surface Wind Speed Based on GNSS-R Forward Model

Abstract

The present invention relates to a method and system for retrieving sea surface wind speed based on the GNSS-R forward model, belonging to the technical field of ocean remote sensing and environmental monitoring. First, by inputting four types of CYGNSS data, a modeled DDM and a Jacobian matrix are obtained. Secondly, feature extraction is carried out. After screening the data quality, multiple features related to the sea surface wind speed are obtained. Then, a spatio-temporal matching with the wind speed is input into the inversion model. Finally, it is concluded that the distributions of the DDM generated by the DDM forward model and the observed DDM are relatively concentrated, and the overall simulation effect is good. The present invention solves the problems that the neural network lacks physical interpretability and the inversion process is complex, and provides new ideas for retrieving soil moisture and other indicators in the future.

Description

Technical Field

[0001] The present invention relates to a method and system for retrieving sea surface wind speed based on a GNSS-R forward model, belonging to the technical field of ocean remote sensing and environmental monitoring. Background Art

[0002] Sea surface wind speed is an important parameter in ocean environmental monitoring, which is of great significance for climate change research, ocean meteorological forecasting, and ocean resource development. Traditional sea surface wind speed measurements usually rely on devices such as buoys or radars, but these methods have problems such as high cost, limited coverage, and complex deployment. GNSS-R is a passive remote sensing technology that retrieves sea surface characteristics by receiving the reflected waves of global navigation satellite signals, and has the advantages of low cost, wide coverage, and all-weather operation. In recent years, significant progress has been made in the application of GNSS-R technology to retrieve sea surface wind speed, especially after combining physical models and data-driven algorithms, the retrieval accuracy has been significantly improved. Therefore, the present invention is proposed to apply GNSS-R technology to sea surface wind speed retrieval. Summary of the Invention

[0003] Aiming at the deficiencies of the prior art, the present invention provides a method and system for retrieving sea surface wind speed based on a GNSS-R forward model, aiming to overcome the problems of high cost and insufficient accuracy in the prior art. Compared with the prior art, the present invention has the advantages of low cost, wide coverage, and strong anti-interference ability. The present invention further improves the accuracy and reliability of wind speed retrieval by introducing the Jacobian matrix for sensitivity analysis and combining the attention mechanism in deep learning.

[0004] Term Explanation:

[0005] CYGNSS: (Cyclone Global Navigation Satellite System) Cyclone Global Navigation Satellite System;

[0006] CYGNSS L1-level GNSS-R data: Cyclone Global Navigation Satellite System level 1 GNSS-R data product;

[0007] DDMI: (Delay Doppler Mapping Instrument) Delay Doppler Imager;

[0008] E2ES: (end-to-end simulator) End-to-end Simulator;

[0009] ERA5: (European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis) European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data.

[0010] The technical solution of the present invention is as follows:

[0011] A method for retrieving sea surface wind speed based on the GNSS-R forward model, the steps are as follows:

[0012] (1) In the GNSS-R receiver, the reflected signal is first cross-correlated with the local copy of the transmitted signal within the delay and Doppler frequency ranges to determine the model input, construct a DDM forward model based on the bistatic scattering model, and calculate the Jacobian matrix;

[0013] (2) Prepare the inversion model data, including CYGNSS L1-level GNSS-R data, ERA5 wind speed and direction, extract the observed variable parameters from the GNSS-R data, perform data quality control, data filtering and model dataset construction, and complete the spatio-temporal matching;

[0014] (3) Combine the Jacobian matrix output by the DDM forward model, construct an inversion model, and evaluate and analyze the inversion results of the model with the introduced Jacobian matrix.

[0015] Preferably according to the present invention, in step (1), the model inputs include four types, namely: geometry, including the positions and velocities of the transmitter, receiver, and specular reflection point;

[0016] Metadata, including sampling time, GPS pseudo-random noise (PRN) code, and specular reflection point index;

[0017] Power, including the effective isotropic radiated power (EIRP) of the GPS transmitter and the antenna parameters of the receiver;

[0018] Wind field, a gridded wind field represented by latitude and longitude coordinates.

[0019] Among them, geometry, metadata, and power all come from CYGNSS level 1 data, and the gridded wind field is constructed from the cross-calibrated multi-platform ocean surface wind field product CCMP dataset.

[0020] Preferably according to the present invention, in step (1), the process of constructing the DDM forward model is:

[0021] The formula of the bistatic scattering model is:

[0022]

[0023] Among them, the integration variables (θ, φ) are two-dimensional position vectors defining the sea surface points, T i represents the coherent integration time of the receiver, λ is the wavelength of the GPS carrier, x(θ, φ) is the wind speed at (θ, φ), P T G T(θ, φ) is the GPS EIRP (effective isotropic radiated power), assumed to be constant over the entire scattering region, G R (θ, φ), G T (θ, φ) represent the gains of the receiver antenna and the transmitter antenna at the surface respectively, P T represents the transmit power, G rain (θ, φ) is the additional rainfall attenuation, represent the path distances from the mean sea surface at the point (θ, φ) to the transmitter and the receiver respectively, and Ω(τ, f, θ, φ) is the surface weighted area for a given delay and Doppler;

[0024] During the simulation process, the semi - empirical Katzberg model is used to convert the wind speed to the surface MSS (Mean Square Slope) sea surface mean square slope. In the E2ES of the CYGNSS's DDMI, Equation (1) is rewritten as a two - dimensional convolution with a delta function to calculate the expected value of the DDM;

[0025] Ω(τ, f, θ, φ) = ∫δ(τ(θ, φ) - τ - τ', f D (θ, φ) - f - f')χ 2 (τ', f')δτ'δf' (2)

[0026] For fast calculation, the fast Fourier transform (FFT) method is applied to the two - dimensional convolution:

[0027] <|Y S (τ, f, x)| 2 >H(τ, f, x)**χ 2 (τ, f) (3)

[0028] In the ambiguity function X 2 (τ, f) and

[0029] H(τ, f, x) = ∫∫h(θ, φ)σ 0 ((θ, φ), x(θ, φ))×δ(τ - τ g (θ, φ))δ(f - f g (τ, φ))dA(θ, φ)

[0030] The δ() between them is the Dirac delta function, and δ() is the brief representation of δ(f - f g (τ, φ)) in the above formula, ** represents the convolution operation, and where h(θ, φ) is the expected value of the power contribution at each position on the ground surface, and H(τ, f, x) is the mapping of this power from the ground surface coordinates to the delay - Doppler coordinates.

[0031] To calculate Equation (3) above, it is necessary to discretize the delay-Doppler space. To model the subsequent Jacobian matrix, the model needs to be inverted together with the observables. Therefore, it is necessary to use the same delay-Doppler coordinates as the observed DDM so that the authenticity of the data can be accurately reflected. The estimated specular point position generated by the DDMI on CYGNSS is very rough and fixed to discrete values. Therefore, the true specular reflection point usually does not lie at the center of a single DDM bin. After use, the fractional accurate specular reflection index and The delay-Doppler coordinates for mobile modeling of DDM are used to align with the observed DDM and the delay-Doppler coordinates of other observables, discretized, and calculated using two-dimensional fast Fourier transform;

[0032] Y[n τ ,n f ,x] = H[n τ ,n f ,x]**Ξ[n τ ,n f (4)

[0033] where Ξ[n τ ,n f = χ 2 (n τ Δτ,n f Δf), n t = 1, 2, …, N τ ,n f = 1, 2, …, N f are indices in the discrete delay-Doppler space, Δτ and Δf are the delay and Doppler increments. For CYGNSS, N τ = 17, N f = 11, Δτ = 0.25T c ,T c is the second period of the GPS C / A code chip, Δf = 1 / (2T i ) = 500Hz, H[n τ ,n f ,x] is the discrete form of the surface integral formula.

[0034] According to the preference of the present invention, in step (1), the Jacobian matrix H is the partial derivative matrix of the forward operator with respect to each wind speed in the input vector x, expressed as follows:

[0035]

[0036] The Jacobian matrix reflects the sensitivity of each sample of the DDM to the wind speed of each surface grid point.

[0037] Preferably according to the present invention, in step (1), in order to improve the physical interpretability of the model and overcome the problem of limited improvement in the accuracy of deep learning methods, by modeling the GNSS-R basic model, a matrix is output that can represent the sensitivity of each sample of the DDM to the wind speed of each surface grid point. Some positions in the matrix are zero because these positions correspond to parts outside the flash area centered on the specular reflection point. The dispersion powers of some functions are very small and are gradually submerged by noise, and the actual values of all Jacobian matrices are zero. High absolute values in the matrix indicate the sensitivity of the DDM to the wind speed at specific positions within the relative delay and Doppler ranges. The two dimensions of the matrix are the DDM values and the surface wind speed, where the DDM is expanded into a vector using delay and Doppler, and the surface wind speed is expanded into a vector using longitude and latitude. To verify the usability of the model, the data of the input variables of the model is obtained from the CYGNSS L1a product at a rate of 1 Hz. The receiver antenna pattern estimated by pre-launch measurement and on-orbit correction is provided by the CYGNSS SOC. At each observation moment, the simulated DDM generated by the forward operator is compared with the DDM obtained from the L1a data. Two metric indicators are used to quantify the consistency between the modeled DDM (the DDM calculated by the DDM forward model) and the observed DDM (the observed DDM is the DDM obtained by processing the CYGNSS L1 data using the delay Doppler imager), which are, respectively, the relative difference of the effective bins and the correlation coefficient between the two DDM effective bins, as follows:

[0038]

[0039] where k is the index of the DDM, y k and h k are the observed DDM and the modeled DDM, i and j respectively represent the indices of the time delay and the Doppler shift, N is the number of effective bins, y k (τ i ,f j ) is the power value of the observed DDM at the time delay τ i and the Doppler shift f j , h k (τ i ,f j ) is the power value of the modeled DDM at the time delay τ i and the Doppler shift f j , τ i is the i-th time delay sampling point, f j is the j-th Doppler shift sampling point;

[0040]

[0041] where ρ k is the correlation coefficient of the k-th DDM, and is the average value of all valid bins for observing DDM and modeling DDM. i and j represent the indices of time delay and Doppler shift respectively, and y k (τ i , f j ) is the power value of the observed DDM at time delay τ i and Doppler shift f j . h k (τ i , f j ) is the power value of the modeled DDM at time delay τ i and Doppler shift f j . τ i is the i-th time delay sampling point, and f j is the j-th Doppler shift sampling point. The more similar the shapes of the observed DDM and the modeled DDM are, the closer the correlation coefficient is to 1.

[0042] Preferably according to the present invention, in step (3), the inversion model architecture is: a CNN module, an auxiliary feature processing module, a feature fusion layer, and an output layer. The CNN module is used to extract the features of a single time delay-Doppler correlation power map. The CNN module consists of 2 convolutional layers, 2 pooling layers, and 1 flattening layer: the 2 convolutional layers each contain 32 3×3 filters and 64 3×3 filters, and are used to perform convolutional operations on the single time delay-Doppler correlation power map to obtain a first feature map; the 2 pooling layers both use max pooling operations and are used to perform max pooling operations on the first feature map to obtain a pooled feature map, with a pooling window size of 2×2 and a stride of 2; the flattening layer is used to flatten the pooled feature map into a first feature sequence with a length of 128;

[0043] The auxiliary feature processing module includes 3 fully connected layers, and the 3 fully connected layers are respectively composed of 64 neurons, 128 neurons, and 256 neurons. The auxiliary feature processing module is used to extract features from auxiliary parameters such as the incident angle, normalized bistatic radar cross-section, and DDM front slope to obtain a second feature sequence;

[0044] The feature fusion layer consists of 2 fully connected layers, and the 2 fully connected layers respectively contain 512 and 256 neurons, and are used to splice and fuse the first feature sequence and the second feature sequence to obtain a fused feature;

[0045] The output layer is a fully connected layer and is used to map the fused feature to a wind speed prediction value;

[0046] The inversion model after introducing the Jacobian matrix is increased by: a sensitivity weighted convolutional layer, a custom attention mechanism, and a dual attention fusion module;

[0047] Finally, define the loss function.

[0048] According to a further preference of the present invention, the sensitivity weighted convolution layer includes a sensitivity mask generation network and a dynamic feature modulation mechanism: the sensitivity mask generation network consists of two layers of convolution networks. The first layer is a 1×1 convolution layer with the number of input channels being the number of channels of the Jacobian matrix and the number of output channels being 32, which is used for channel mapping. The second layer is a 3×3 convolution layer with the number of input channels being 32 and the number of output channels being the same as the number of original feature channels, which is used to capture spatial correlation. The dynamic feature modulation mechanism includes a learnable scaling factor for dynamically adjusting the generated sensitivity mask and then weighting the input features.

[0049] Generate a sensitivity mask from the Jacobian matrix through two layers of convolution networks:

[0050] FM(J) = Conv2d(ReLU(J, C1, 1x1), C2, kxk) (8)

[0051] Where J is the input Jacobian matrix, C1 is the number of channels of the first convolution layer, set to 32, C2 is the number of channels of the second convolution layer, set to the number of input channels, k is the size of the second convolution kernel, set to 3. ReLU is the activation function.

[0052] Dynamic feature modulation mechanism:

[0053] X' = X * (1 + α * FM(J)) (9)

[0054] Where X is the input feature, α is the learnable scaling factor, initialized to 1.0, and * represents element-wise multiplication.

[0055] According to a preference of the present invention, the custom attention mechanism consists of a multi-head attention structure and Jacobian-guided attention weights: the multi-head attention structure divides the input features into 4 attention heads, and each attention head independently calculates the attention weights to capture the feature correlations in different subspaces. The Jacobian-guided attention weights convert the Jacobian matrix information into attention modulation coefficients through a two-layer perceptron. The output dimension of the first layer is 64, and the output dimension of the second layer is the number of attention heads, which is used to dynamically adjust the attention distribution.

[0056] Jacobian row-guided attention weights:

[0057] Att_weight = σ(MLP(J)) (10)

[0058] Att_out = Attention(Q, K, V) * (1 + β * Att_weight) (11)

[0059] Among them, σ is the sigmoid activation function, β is a scientific parameter, Q (Query), K (Key), and V (Value) are three key components in the attention mechanism, which are the query matrix, the key matrix, and the value matrix respectively. MLP is a multi-layer perceptron: the number of parameters in the first layer = input_dim × 64 + 64 (bias term), the number of parameters in the second layer = 64 × num_heads + num_heads (bias term), input_dim is the feature dimension of the Jacobian matrix, 64 neurons are selected in the middle for feature transformation and dimensionality increase, and the output is the attention weight of the num_heads dimension. The main role of MLP is to convert the information of the Jacobian matrix into the weight coefficients of each head in the attention mechanism, so as to achieve attention modulation based on sensitivity.

[0060] Preferably according to the present invention, the dual attention fusion module includes channel attention and spatial attention: the channel attention calculates the importance weight of each channel through global average pooling and a two-layer perceptron. The first-layer perceptron reduces the number of channels by 16 times, and the second-layer perceptron restores the original number of channels;

[0061] The spatial attention uses a 7×7 convolutional layer to generate an attention map in the spatial dimension, and the number of output channels is 1. The outputs of the two attentions act on the original features through element-wise multiplication to achieve feature enhancement.

[0062] Preferably according to the present invention, the loss function consists of a basic loss term, a sensitivity weighting term, and a regularization term: the basic loss term uses the mean square error; the sensitivity weighting term calculates the sample weight based on the absolute value of the Jacobian matrix; the regularization term uses L2 regularization to prevent overfitting;

[0063] The loss function is:

[0064]

[0065] Among them, Sensitivity_Loss is the sensitivity weighted loss, |J| is the absolute value of the Jacobian matrix, L2_reg is the L2 regularization term, ||W||_2 is the L2 norm of the weight, that is, the square root of the sum of the squares of all weights, W is the weight parameter of the model, λ1 is the weight coefficient of the sensitivity loss, set to 0.1, λ2 is the weight coefficient of L2 regularization, set to 0.001, and MSE is the mean square error:

[0066] MSE = 1 / n * ∑(y_pred - y_true) 2 (13)

[0067] Among them, y_pred is the wind speed value predicted by the model, y_true is the true wind speed value, and n is the number of samples.

[0068] Preferably according to the present invention, in step (3), the model inversion result is quantified by the following formula:

[0069]

[0070] where is the predicted wind speed, v is the reference wind speed from ERA5, the overall RMSE and bias are determined based on a total of n samples of all models, the root mean square error is applicable to all inversion architectures and the MVE traditional model, RMSE (root mean square error) reflects the average deviation between the model prediction value and the true value, the smaller the value, the higher the model prediction accuracy, Bias (bias) reflects the systematic deviation of the prediction value relative to the true value, and is used to judge whether the model has a tendency of overestimation or underestimation. RMSE reflects the overall accuracy of the prediction, and Bias reflects the systematic error of the prediction.

[0071] A sea surface wind speed inversion system based on the GNSS-R forward model includes:

[0072] A modeling module: used to perform cross-correlation between the reflected signal and the local copy of the transmitted signal within the delay and Doppler frequency range in the GNSS-R receiver, determine the model input, construct a DDM forward model based on the bistatic scattering model, and calculate the Jacobian matrix;

[0073] A data preparation module, used to prepare the inversion model data, including CYGNSS L1-level GNSS-R data, ERA5 wind speed and direction, extract the observed variable parameters from the GNSS-R data, perform data quality control, data filtering, and model dataset construction, and complete the spatio-temporal matching;

[0074] An inversion construction module, used to combine the DDM forward model output Jacobian matrix, construct an inversion model, and evaluate and analyze the model inversion result introduced with the Jacobian matrix.

[0075] The beneficial effects of the present invention are as follows:

[0076] First, by inputting four types of CYGNSS data, the present invention obtains a modeled DDM (Delay-Doppler Map) and a Jacobian matrix. Secondly, feature extraction is carried out. After screening the data quality, multiple features related to the sea surface wind speed are obtained. Then, spatio-temporal matching with the wind speed is performed and input into the inversion model. Finally, it is concluded that the distributions of the DDM generated by the DDM forward model and the observed DDM are relatively concentrated, and the overall simulation effect is good. Moreover, the relative error values of the samples basically remain below 20%. The distribution of the correlation coefficients is mainly concentrated between 0.94 and 0.98, and most values are close to 1. Only a small number of sample values are below 0.92, and the correlation coefficient will decrease in the case of high wind speeds. In summary, it can be concluded that the inversion model performs well, and the performance of the model on most samples is very close to the observed data. The modeled DDM is basically consistent with the DDM measured by CYGNSS. The RMSE of the inversion model has also increased by 16% compared with the original model without the introduction of the Jacobian matrix, from 1.875 m / s to 1.572 m / s. The present invention solves the problems that the neural network lacks physical interpretability and the inversion process is complex, and provides new ideas for the subsequent inversion of soil moisture and other indicators. BRIEF DESCRIPTION OF THE DRAWINGS

[0077] Figure 1 is a schematic flow chart of the present invention;

[0078] Figure 2 is a comparison diagram of the modeled DDM and the observed DDM of the present invention, wherein, Figure 2 in (a) is the observed DDM diagram, Figure 2 in (b) is the modeled DDM diagram. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0079] The present invention will be further described below by way of examples in conjunction with the drawings, but is not limited thereto.

[0080] Example 1:

[0081] As Figure 1 shown, this embodiment provides a method for inverting the sea surface wind speed based on the GNSS-R forward model, and the steps are as follows:

[0082] (1) In the GNSS-R receiver, the reflected signal is first cross-correlated with the local copy of the transmitted signal within the delay and Doppler frequency ranges to determine the model input, construct a DDM forward model based on the bistatic scattering model, and calculate the obtained Jacobian matrix;

[0083] The model input includes four types, namely: geometry, including the positions and velocities of the transmitter, receiver, and specular reflection point;

[0084] Metadata, including sampling time, GPS pseudo-random noise (PRN) code, and specular reflection point index;

[0085] Power, including the effective isotropic radiated power (EIRP) of the GPS transmitter and the antenna parameters of the receiver;

[0086] Wind field, a gridded wind field represented by latitude and longitude coordinates.

[0087] Among them, geometry, metadata, and power all come from CYGNSS level 1 data, and the gridded wind field is constructed from the Cross-Calibrated Multi-Platform Ocean Surface Wind Fields Product CCMP dataset.

[0088] The process of constructing the DDM forward model is as follows:

[0089] The bistatic scattering model formula is:

[0090]

[0091] Among them, the integration variables (θ, φ) are the two-dimensional position vectors defining the sea surface points, T i represents the coherent integration time of the receiver, λ is the wavelength of the GPS carrier, x(θ, φ) is the wind speed at (θ, φ), P T G T (θ, φ) is the effective isotropic radiated power (EIRP) of the GPS, assumed to be constant throughout the scattering region, G R (θ, φ), G T (θ, φ) represent the gains of the receiver antenna and the transmitter antenna at the surface respectively, P T represents the transmit power, G rain (θ, φ) is the additional rainfall attenuation, represent the path distances from the mean sea surface at point (θ, φ) to the transmitter and the receiver respectively, and Ω(τ, f, θ, φ) is the surface weighted area for a given delay and Doppler;

[0092] During the simulation process, the semi-empirical Katzberg model is used to convert the wind speed to the surface Mean Square Slope (MSS). In the E2ES of the DDMI in CYGNSS, equation (1) is rewritten as a two-dimensional convolution with a delta function to calculate the expected value of the DDM;

[0093] Ω(τ, f, θ, φ) = ∫δ(τ(θ, φ) - τ - τ', f D (θ, φ) - f - f')χ 2 (τ', f')δτ'δf' (2)

[0094] For fast calculation, the fast Fourier transform (FFT) method is applied to two-dimensional convolution:

[0095] <|Y S (τ,f,x)| 2 >=H(τ,f,x)**χ 2 (τ,f) (3)

[0096] In the ambiguity function X 2 (τ,f) and

[0097] H(τ,f,x)=∫∫h(θ,φ)σ 0 ((θ,φ),x(θ,φ))×δ(τ-τ g (θ,φ))δ(f-f g (τ,φ))dA(θ,φ)

[0098] The δ() between them is the Dirac delta function, and δ() is the brief representation of δ(f-f g (τ,φ)) in the above formula. ** represents the convolution operation, and where h(θ,φ) is the expected value of the power contribution at each position on the ground surface, and H(τ,f,x) is the mapping of this power from the ground surface coordinates to the delay-Doppler coordinates.

[0099] To calculate equation (3) above, it is necessary to discretize the delay-Doppler space. To model the subsequent Jacobian matrix, the model inversion needs to be carried out together with the observables. Therefore, it is necessary to use the same delay-Doppler coordinates as the observed DDM so that the authenticity of the data can be accurately reflected. The mirror point position estimation generated by the DDMI on CYGNSS is very rough and fixed as discrete values. Therefore, the true mirror reflection point usually does not lie at the center of a single DDM bin. The fractional accurate mirror reflection index generated later and are used to move and model the delay-Doppler coordinates of the DDM to align with the delay-Doppler coordinates of the observed DDM and other observables, perform discretization processing, and use two-dimensional fast Fourier transform for calculation;

[0100] Y[n τ ,n f ,x]=H[n τ ,n f ,x]**Ξ[n τ ,n f (4)

[0101] where Ξ[n τ ,n f =χ 2 (n τ Δτ,nf Δf), n t = 1, 2, …, N τ , n f = 1, 2, …, N f are indices in the discrete delay - Doppler space, Δτ and Δf are delay and Doppler increments. For CYGNSS, N τ = 17, N f = 11, Δτ = 0.25T c , T c is the second - period of the GPS C / A code chip, Δf = 1 / (2T i ) = 500Hz, H[n τ , n f , x] is the discrete form of the surface integral formula.

[0102] The Jacobian matrix H is the matrix of partial derivatives of the forward operator with respect to each wind speed in the input vector x and is expressed as follows:

[0103]

[0104] The Jacobian matrix reflects the sensitivity of each sample of the DDM to the wind speed at each surface grid point.

[0105] To improve the physical interpretability of the model and overcome the problem of limited accuracy improvement in deep - learning methods, by modeling the GNSS - R basic model, the output can be a matrix representing the sensitivity of each sample of the DDM to the wind speed at each surface grid point. Some positions in the matrix are zero because these positions correspond to parts outside the flash region centered on the specular reflection point, and the dispersion powers of some functions are very small and gradually submerged by noise. The actual values of all Jacobian matrices are zero. High absolute values in the matrix indicate the sensitivity of the DDM to the wind speed at a specific position within the relative delay and Doppler ranges. The two dimensions of the matrix are the DDM values and the surface wind speeds, where the DDM is expanded into a vector using delay and Doppler, and the surface wind speeds are expanded into a vector using longitude and latitude. To verify the usability of the model, the data of the input variables of the model are obtained from CYGNSS L1a products at a rate of 1Hz. The receiver antenna pattern estimated through pre - launch measurement and on - orbit correction is provided by the CYGNSS SOC. At each observation moment, the simulated DDM generated by the forward operator is compared with the DDM obtained from the L1a data. Two metric indices are used to quantify the consistency between the modeled DDM (the DDM calculated through the DDM forward model) and the observed DDM (the observed DDM is the DDM obtained by processing the CYGNSS L1 data using the delay - Doppler imager), which are, respectively, the relative difference of the valid bins and the correlation coefficient between the two DDM valid bins, as follows:

[0106]

[0107] where k is the index of the DDM, y k and h k are the observed DDM and the modeled DDM, i and j represent the indices of time delay and Doppler shift respectively, N is the number of effective bins, y k (τ i , f j ) is the power value of the observed DDM at time delay τ i and Doppler shift f j , h k (τ i , f j ) is the power value of the modeled DDM at time delay τ i and Doppler shift f j , τ i is the i-th time delay sampling point, f j is the j-th Doppler shift sampling point;

[0108]

[0109] where ρ k is the correlation coefficient of the k-th DDM, and are the averages of all effective bins of the observed DDM and the modeled DDM, i and j represent the indices of time delay and Doppler shift respectively, y k (τ i , f j ) is the power value of the observed DDM at time delay τ i and Doppler shift f j , h k (τ i , f j ) is the power value of the modeled DDM at time delay τ i and Doppler shift f j , τ i is the i-th time delay sampling point, f j is the j-th Doppler shift sampling point. The more similar the shapes of the observed DDM and the modeled DDM are, the closer the correlation coefficient is to 1.

[0110] The comparison results are as Figure 2As shown, the simulated DDM and Jacobian matrix are obtained through the DDM forward model. Evaluating the model reveals that the relative difference of the DDM is less than 30%, and the correlation coefficient is greater than 90%, ensuring the consistency of the model and also demonstrating the corresponding relationship between the Jacobian matrix and the DDM under different wind speeds. Under ideal circumstances, the main reason for the result difference in the DDM forward model is the difference between the model wind field and the real wind. Additionally, there are some other uncertain factors that can affect the performance of the model, including: at low wind speeds in medium and high wind speed conditions, local winds play a dominant role, so the influence of swells can be ignored. However, at low wind speeds, the influence of swells is obvious. For the problem of inaccurate specular reflection point positions, more accurate specular point calculation methods and a higher-fidelity sea surface model can be used to calculate the specular reflection point delay and Doppler. For the inaccurately estimated GPS transmitter EIRP, its estimation accuracy is affected by transmitter power changes and zenith antenna pattern calibration. Subsequently, the accuracy of the forward model can be improved by perfecting these aspects.

[0111] The modeled DDM and Jacobian matrix generated by the DDM forward model, and the specific variable names in the CYGNSS level 1 data involved are as follows in the table:

[0112]

[0113] Modeled according to the GPS bistatic radar scattering model, the DDM (in units of the absolute power at the receiving end) is related to the sea surface wind field.

[0114] The ground truth wind speed uses the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis data. The time resolution of the ERA5 wind speed product is 1 h, and the grid resolution is 0.5×0.5. The longitude of the ERA5 data is converted to the range of 0 - 360. An interpolation function is created for the wind field file, and then interpolation is performed in the ERA5 data using the longitude and latitude of CYGNSS. Finally, the obtained wind speed is saved to a new variable in the CYGNSS data file. Through time matching and spatial interpolation in the whole process, it can be ensured that each observation point of CYGNSS can obtain the corresponding ERA5 wind speed. High-quality DDM observation data is essential for the model to correctly and intuitively understand how to map different ocean wind speeds to their corresponding observations. Therefore, it is essential to eliminate the low-quality measurement values in the data, specifically as follows: the uncertainty of the uncertainty of the bistatic radar cross-section (ddm_brcs_uncert) of BRCS is less than 1;

[0115] The attitude status of the nanosatellite tracker is "OK", that is

[0116] nst_att_status = 0, the attitude status of the navigation science terminal;

[0117] The absolute value of the spacecraft roll is between 1 degree and 30 degrees, the spacing is between 1 and 10, or the yaw angle is between 1 and 5, indicated by the quality_flags quality identification of each DDM

[0118] For the specular reflection point (sp_rx_gain), the receiving antenna gain in the direction of the receiving antenna gain of the specular reflection point is greater than 0 dBi.

[0119] For DDM (prn_fig_of_merit), the figure of merit (FOM) of the range correction gain (RCG) for PRN selection quality assessment is greater than 0

[0120] For LES (ddm_les), the waveform front slope of the DDM is greater than 0

[0121] For Zenith (direct), the signal-to-noise ratio (direct_signal_snr) of the zenith (direct) signal is greater than 0 dB

[0122] In addition, since CYGNSS measurements or ERA5 estimates may also be contaminated by unexpected errors and these errors usually do not affect the main group, samples outside the 95% confidence interval are excluded.

[0123] (2) Prepare the inversion model data, including CYGNSS L1-level GNSS-R data, ERA5 wind speed and direction, extract the observed variable parameters from the GNSS-R data, perform data quality control, data filtering, and construct the model dataset to complete the spatio-temporal matching;

[0124] (3) Combine the DDM forward model to output the Jacobian matrix, construct the inversion model, and evaluate and analyze the inversion results of the model with the introduced Jacobian matrix.

[0125] The inversion model architecture is: a CNN module, an auxiliary feature processing module, a feature fusion layer, and an output layer. The CNN module is used to extract the features of a single delay-Doppler correlation power map. The CNN module consists of 2 convolutional layers, 2 pooling layers, and 1 flattening layer: the 2 convolutional layers each contain 32 3×3 filters and 64 3×3 filters, which are used to perform convolution operations on the single delay-Doppler correlation power map to obtain the first feature map; both 2 pooling layers use max pooling operations to perform max pooling operations on the first feature map to obtain the pooled feature map, with a pooling window size of 2×2 and a stride of 2; the flattening layer is used to flatten the pooled feature map into a first feature sequence with a length of 128;

[0126] The auxiliary feature processing module contains 3 fully connected layers, and the 3 fully connected layers are respectively composed of 64 neurons, 128 neurons, and 256 neurons. The auxiliary feature processing module is used to extract features from auxiliary parameters such as the incident angle, normalized bistatic radar cross section, and DDM front slope to obtain a second feature sequence;

[0127] The feature fusion layer consists of 2 fully connected layers. The 2 fully connected layers respectively contain 512 and 256 neurons, and are used to splice and fuse the first feature sequence and the second feature sequence to obtain a fused feature;

[0128] The output layer is a fully connected layer, which is used to map the fused feature to a wind speed prediction value;

[0129] The inversion model after introducing the Jacobian matrix is increased with: a sensitivity weighted convolution layer, a custom attention mechanism, and a dual attention fusion module;

[0130] Finally, a loss function is defined.

[0131] The sensitivity weighted convolution layer contains a sensitivity mask generation network and a dynamic feature modulation mechanism: The sensitivity mask generation network consists of two layers of convolutional networks. The first layer is a 1×1 convolutional layer, the number of input channels is the number of channels of the Jacobian matrix, and the number of output channels is 32, which is used for channel mapping; the second layer is a 3×3 convolutional layer, the number of input channels is 32, and the number of output channels is the same as the number of original feature channels, which is used to capture spatial correlation; the dynamic feature modulation mechanism contains a learnable scaling factor, which is used to dynamically adjust the generated sensitivity mask, and then weight the input features;

[0132] Generate a sensitivity mask from the Jacobian matrix through two layers of convolutional networks:

[0133] FM(J) = Conv2d(ReLU(J, Cl, 1x1), C2, kxk) (8)

[0134] Where J is the input Jacobian matrix, C1 is the number of channels of the first convolutional layer, set to 32, C2 is the number of channels of the second convolutional layer, set to the number of input channels, k is the size of the second convolutional kernel, set to 3. ReLU is the activation function;

[0135] Dynamic feature modulation mechanism:

[0136] X' = X * (1 + α * FM(J)) (9)

[0137] Where X is the input feature, α is a learnable scaling factor, initialized to 1.0, and * represents element-wise multiplication.

[0138] Preferably according to the present invention, the custom attention mechanism consists of a multi-head attention structure and Jacobian-guided attention weights: the multi-head attention structure divides the input features into 4 attention heads, and each attention head independently calculates the attention weights to capture the feature correlations in different subspaces; the Jacobian-guided attention weights convert the Jacobian matrix information into attention modulation coefficients through a two-layer perceptron. The output dimension of the first layer is 64, and the output dimension of the second layer is the number of attention heads, which is used to dynamically adjust the attention distribution;

[0139] Jacobian row-guided attention weights:

[0140] Att_weight = σ(MLP(J)) (10)

[0141] Att_out = Attention(Q, K, V) * (1 + β * Att_weight) (11)

[0142] Where, σ is the sigmoid activation function, β is a scientific parameter, Q (Query), K (Key), and V (Value) are three key components in the attention mechanism, which are the query matrix, the key matrix, and the value matrix respectively. MLP is a multi-layer perceptron: the number of parameters in the first layer = input_dim × 64 + 64 (bias term), the number of parameters in the second layer = 64 × num_heads + num_heads (bias term), input_dim is the feature dimension of the Jacobian matrix, 64 neurons are selected in the middle for feature transformation and dimension elevation, and the output is the attention weights of the num_heads dimension. The main role of MLP is to convert the information of the Jacobian matrix into the weight coefficients of each head in the attention mechanism, so as to achieve sensitivity-based attention modulation.

[0143] The dual attention fusion module includes channel attention and spatial attention: channel attention calculates the importance weights of each channel through global average pooling and a two-layer perceptron. The first-layer perceptron reduces the number of channels by 16 times, and the second-layer perceptron restores the original number of channels;

[0144] Spatial attention uses a 7×7 convolutional layer to generate an attention map in the spatial dimension, and the output number of channels is 1. The outputs of the two kinds of attention act on the original features through element-wise multiplication to achieve feature enhancement.

[0145] The loss function consists of a basic loss term, a sensitivity weighting term, and a regularization term: the basic loss term uses mean squared error; the sensitivity weighting term calculates the sample weights based on the absolute value of the Jacobian matrix; the regularization term uses L2 regularization to prevent overfitting;

[0146] The loss function is:

[0147]

[0148] Among them, Sensitivity_Loss is the sensitivity weighted loss, |J| is the absolute value of the Jacobian matrix, L2_reg is the L2 regularization term, ||W||_2 is the L2 norm of the weights, which is the square root of the sum of the squares of all weights, W is the weight parameter of the model, λ1 is the weight coefficient of the sensitivity loss, set to 0.1, λ2 is the weight coefficient of the L2 regularization, set to 0.001, and MSE is the mean squared error:

[0149] MSE = 1 / n * Σ(y_pred - y_true) 2 (13)

[0150] Among them, y_pred is the wind speed value predicted by the model, y_true is the true wind speed value, and n is the number of samples.

[0151] The evaluation of the inversion model was carried out during the test period from May 1 to July 2017. Since the research focus is on the effect of introducing physical constraints, several one-dimensional variable models of DDM variables were selected for effect feedback. The inversion results of the model were quantified by the following formula:

[0152]

[0153] Among them, is the predicted wind speed, v is the reference wind speed from ERA5. The overall RMSE and bias are determined based on a total of n samples of all models. The root mean square error is applicable to all inversion architectures and the MVE traditional model. RMSE (root mean square error) reflects the average deviation between the model prediction value and the true value. The smaller the value, the higher the prediction accuracy of the model. Bias (deviation) reflects the systematic deviation of the prediction value relative to the true value and is used to judge whether there is a trend of overestimation or underestimation in the model. RMSE reflects the overall accuracy of the prediction, and Bias reflects the systematic error of the prediction.

[0154] The research data mainly comes from the CYGNSS project and is verified using the ERA5 wind speed product. The accuracy changes of the model inversion before and after introducing the Jacobian matrix were mainly compared. The results show that the inverted wind of the new model shows good consistency with the ERA5 wind speed product.

[0155] Example 2:

[0156] This example provides a sea surface wind speed inversion system based on the GNSS-R forward model, including:

[0157] Modeling module: used in the GNSS-R receiver to cross-correlate the reflected signal with the local copy of the transmitted signal within the range of delay and Doppler frequency, determine the model input, construct the DDM forward model based on the bistatic scattering model, and calculate the Jacobian matrix;

[0158] Data preparation module: used to prepare the inversion model data, including CYGNSS L1-level GNSS-R data, ERA5 wind speed and direction, extract the observed variable parameters from the GNSS-R data, perform data quality control, data filtering and model dataset construction, and complete the spatio-temporal matching;

[0159] Inversion construction module: used to combine the Jacobian matrix output by the DDM forward model, construct the inversion model, and evaluate and analyze the inversion results of the model with the introduced Jacobian matrix.

[0160] The description of the above embodiments is only used to help understand the method of the present invention and its core idea; at the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manner and application scope. In summary, the content of this specification should not be construed as a limitation to the present invention.

[0161] The above embodiments are used to explain the present invention, rather than limiting the present invention. Any modifications and changes made to the present invention within the spirit and scope of the claims of the present invention fall within the protection scope of the present invention.

Claims

1. The sea surface wind speed inversion method based on the GNSS-R forward model is characterized by: Here are the steps: (1) In the GNSS-R receiver, the reflected signal is first cross-correlated with a local replica of the transmitted signal within the delay and Doppler frequency range to determine the model input, build a DDM forward model based on the bistatic scattering model, and calculate the Jacobian matrix; (2) Prepare inversion model data, including CYGNSS L1 GNSS-R data, ERA5 wind speed and direction, extract observation variable parameters from GNSS-R data, perform data quality control, data filtering, and model data set construction to complete spatiotemporal matching; (3) Combine the Jacobian matrix output by the DDM forward model to construct an inversion model, and evaluate and analyze the inversion results of the model with the Jacobian matrix introduced; The inversion model architecture is: CNN module, auxiliary feature processing module, feature fusion layer and output layer. The CNN module is used to extract the features of a single delay-Doppler correlation power map to obtain the first feature sequence. The CNN module consists of 2 convolutional layers, 2 pooling layers and 1 flattening layer. The auxiliary feature processing module includes three fully connected layers, and is used to extract features from auxiliary parameters to obtain a second feature sequence; The feature fusion layer is composed of two fully connected layers, which are used to concatenate and fuse the first feature sequence and the second feature sequence to obtain a fusion feature; The output layer is a fully connected layer, which is used to map the fused features into wind speed prediction values; The inversion model after the introduction of the Jacobian matrix adds: sensitivity weighted convolution layer, custom attention mechanism and dual attention fusion module; Finally, define the loss function; The sensitivity weighted convolution layer includes a sensitivity mask generation network and a dynamic feature modulation mechanism: the sensitivity mask generation network consists of two layers of convolutional networks. The first layer is a 1×1 convolutional layer, the number of input channels is the number of Jacobian matrix channels, and the number of output channels is 32, which is used for channel mapping; the second layer is a 3×3 convolutional layer, the number of input channels is 32, and the number of output channels is the same as the number of original feature channels, which is used to capture spatial correlation; The dynamic feature modulation mechanism includes a scaling factor to dynamically adjust the generated sensitivity mask to weight the input features; The sensitivity mask is generated from the Jacobian matrix through a two-layer convolutional network: FM(J)=Conv2d(ReLU(J,Cl,1x1),C2,k×k) (8) Where J is the input Jacobian matrix, C1 is the number of convolution channels in the first layer, which is set to 32, C2 is the number of convolution channels in the second layer, which is set to the number of input channels, k is the size of the convolution kernel in the second layer, which is set to 3, and ReLU is the activation function; Dynamic feature modulation mechanism: X'=X*(1+α*FM(J)) (9) Where X is the input feature, α is a learnable scaling factor, initialized to 1.0, and * represents term-by-term multiplication; The custom attention mechanism consists of a multi-head attention structure and Jacobi-guided attention weights: the multi-head attention structure divides the input features into four attention heads, and each attention head independently calculates the attention weight to capture the feature correlation of different subspaces; the Jacobi-guided attention weights convert the Jacobi matrix information into attention modulation coefficients through a two-layer perceptron. The output dimension of the first layer is 64, and the output dimension of the second layer is the number of attention heads, which is used to dynamically adjust the attention distribution; Jacobian-guided attention weights: Att_weight = σ(MLP(J)) (10) Att_out=Attention(Q,K,V)*(1+β*Att_weight) (11) Among them, σ is the sigmoid activation function, β is the scientific parameter, Q, K, and V are the three key components in the attention mechanism, namely the query matrix, key matrix, and value matrix, respectively, and MLP is a multi-layer perceptron; The dual attention fusion module includes channel attention and spatial attention: channel attention calculates the importance weight of each channel through global average pooling and two-layer perceptron. The first layer of perceptron reduces the number of channels by 16 times, and the second layer of perceptron restores the original number of channels. Spatial attention uses a 7×7 convolutional layer to generate an attention map of the spatial dimension, with an output channel of 1. The outputs of the two attentions are applied to the original features through element-by-element multiplication to achieve feature enhancement.

2. The sea surface wind speed inversion method based on the GNSS-R forward model according to claim 1, characterized in that: In step (1), the model inputs include four types, namely: geometry, including the positions and velocities of the transmitter, receiver, and specular reflection points; Metadata, including sampling time, GPS pseudo-random noise code, and specular reflection point index; Power, including the effective isotropic radiated power of the GPS transmitter and the antenna parameters of the receiver; Wind field, a gridded wind field expressed in latitude and longitude coordinates.

3. The sea surface wind speed inversion method based on the GNSS-R forward model as claimed in claim 2, characterized in that: In step (1), the DDM forward model construction process is: The bistatic scattering model formula is: The integral variables (θ, φ) are the two-dimensional position vectors defining the sea surface point, T i represents the coherent integration time of the receiver, λ is the wavelength of the GPS carrier, x(θ,φ) is the wind speed at (θ,φ), P T G T (θ, φ) is the GPS EIRP equivalent isotropic radiated power, which is assumed to remain constant throughout the scattering region. R (θ,φ),G T (θ, φ) represent the gain of the receiver antenna and the transmitter antenna on the surface, respectively, P T Indicates the transmit power, G rain (θ,φ) is the additional rain attenuation, denote the path distances from the average sea surface at point (θ, φ) to the transmitter and receiver, respectively, and Ω(τ, f, θ, φ) is the surface weighted area for a given delay and Doppler; During the simulation, the semi-empirical Katzberg model is used to convert the wind speed into the surface MSS. In E2ES in the DDMI of CYGNSS, equation (1) is rewritten as a two-dimensional convolution with a delta function to calculate the expected value of the DDM; Ω(τ,f,θ,φ)=∫δ(τ(θ,φ)-τ-τ',f D (θ,φ)-f-f')χ 2 (τ',f')δτ'δf' (2) For fast computation, the Fast Fourier Transform (FFT) method is applied to the two-dimensional convolution: <|Y S (τ,f,x)| 2 >=H(τ,f,x)**x 2 (τ,f) (3) In the ambiguity function χ 2 (τ,f) and H(τ,f,x)=∫∫h(θ,φ)σ 0 ((θ,φ),x(θ,φ))×δ(τ-τ g (θ,φ))δ(ff g (τ,φ))dA(θ,φ) The δ() between is the Dirac delta function, ** represents the convolution operation, and Where h(θ,φ) is the expected value of the power contribution at each location on the surface, and H(τ,f,x) is the mapping of this power from surface coordinates to delay-Doppler coordinates; To calculate the above equation (3), the delay-Doppler space needs to be discretized and calculated using a two-dimensional fast Fourier transform; Y[n τ ,n f ,x]=H[n τ ,n f ,x]**Ξ[n τ ,n f ] (4) Among them, Ξ[n τ ,n f ]=χ 2 (n τ Δτ,n f Δf), n t =1,2,…,N τ , n f =1,2,…,N f is the index in discrete delay-Doppler space, Δτ and Δf are the delay and Doppler increments, and for CYGNSS, N τ =17,N f =11, Δτ=0.25T c , T c is the second period of GPSC / A code chip, Δf=1 / (2T i )=500Hz,H[n τ ,n f ,x] is the discrete form of the surface integral.

4. The sea surface wind speed inversion method based on the GNSS-R forward model as claimed in claim 3, characterized in that: In step (1), the Jacobian matrix H is the partial derivative matrix of the forward operator with respect to each wind speed in the input vector x, expressed as follows: The Jacobian matrix reflects the sensitivity of each sample of the DDM to the wind speed at each surface grid point.

5. The sea surface wind speed inversion method based on the GNSS-R forward model as claimed in claim 4, characterized in that: In step (1), two metrics are used to quantify the consistency between the modeled DDM and the observed DDM, namely, the relative difference of the valid bins and the correlation coefficient between the valid bins of the two DDMs, as follows: Where k is the index of DDM, y k and h k are the observed DDM and the modeled DDM, i and j represent the indexes of delay and Doppler shift, respectively, N is the number of valid bins, and y k (τ i ,f j ) is the observed DDM at a delay τ i and Doppler frequency shift f j The power value at k (τ i ,f j ) is the modeling DDM at the time delay τ i and Doppler frequency shift f j The power value at τ i is the i-th delayed sampling point, f j is the jth Doppler frequency shift sampling point; Among them, ρ k is the correlation coefficient of the kth DDM, and is the average of all valid bins of the observed DDM and the modeled DDM, i and j represent the indexes of delay and Doppler shift, respectively, and y k (τ i ,f j ) is the observed DDM at a delay τ i and Doppler frequency shift f j The power value at k (τ i ,f j ) is the modeling DDM at the time delay τ i and Doppler frequency shift f j The power value at τ i is the i-th delayed sampling point, f j is the jth Doppler shift sampling point. The more similar the shapes of the observed DDM and the modeled DDM are, the closer the correlation coefficient is to 1.

6. The sea surface wind speed inversion method based on the GNSS-R forward model as claimed in claim 5, characterized in that: The loss function consists of a basic loss term, a sensitivity weighted term, and a regularization term: the basic loss term uses mean square error; the sensitivity weighted term calculates the sample weight based on the absolute value of the Jacobian matrix; the regularization term uses L2 regularization to prevent overfitting; The loss function is: Among them, Sensitivity_Loss is the sensitivity weighted loss, |J| is the absolute value of the Jacobian matrix, L2_reg is the L2 regularization term, ||W||_2 is the L2 norm of the weight, that is, the square root of the sum of all squares, W is the weight parameter of the model, λ1 is the weight coefficient of the sensitivity loss, set to 0.1, λ2 is the weight coefficient of the L2 regularization, set to 0.001, and MSE is the mean square error: MSE=1 / n*∑(y_pred-y_true) 2 (13) Among them, y_pred is the wind speed value predicted by the model, y_true is the true wind speed value, and n is the number of samples.

7. The sea surface wind speed inversion method based on the GNSS-R forward model as claimed in claim 6, characterized in that: In step (3), the model inversion results are quantified by the following formula: in, is the predicted wind speed, v is the reference wind speed from ERA5, RMSE reflects the average deviation between the model prediction value and the true value. The smaller the value, the higher the model prediction accuracy. Bias reflects the systematic deviation of the predicted value relative to the true value, which is used to determine whether the model has a tendency to overestimate or underestimate. RMSE reflects the overall accuracy of the prediction, and Bias reflects the systematic error of the prediction.

8. A sea surface wind speed inversion system based on a GNSS-R forward model, applied to the sea surface wind speed inversion method based on a GNSS-R forward model as claimed in claim 1, characterized in that: include: Modeling module: used in the GNSS-R receiver to cross-correlate the reflected signal with the local copy of the transmitted signal within the delay and Doppler frequency range, determine the model input, construct the DDM forward model based on the bistatic scattering model, and calculate the Jacobian matrix; Data preparation module, used to prepare inversion model data, including CYGNSS L1 GNSS-R data, ERA5 wind speed and direction, extract observation variable parameters from GNSS-R data, perform data quality control, data filtering and model data set construction, and complete time-space matching; The inversion construction module is used to combine the Jacobian matrix output by the DDM forward model, construct the inversion model, and evaluate and analyze the inversion results of the model with the Jacobian matrix introduced.

Images