Joint retrieval method of ocean swell height and wind wave height from space-borne GNSS-R multi-modal data fusion

Abstract

The application discloses a joint inversion method for ocean swell height and wind wave height based on satellite-borne GNSS-R multi-modal data fusion, belongs to the cross-research field of GNSS reflection measurement and deep learning, and comprises the following steps: S1, obtaining original data; S2, processing the original data; S3, dividing the data set processed in the step S2 into training data, verification data and test data; S4, constructing a joint inversion model for ocean swell height and wind wave height; S5, adaptively weighting the two targets of ocean swell height and wind wave height; S6, inputting the test data into the joint inversion model for testing; and S7, evaluating the performance of the model.The model has the beneficial effects that the model can directly output the two targets at one time, and at the same time, the inversion accuracy of the two targets is ensured to a certain extent, the calculation efficiency is greatly saved, and the inversion effect of the model is not inferior to that of a single output model, and the economic benefit is achieved.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of GNSS reflection measurement and deep learning cross-research, in particular to a joint inversion method for ocean swell height and wind wave height based on multi-modal data fusion of spaceborne GNSS-R. BACKGROUND

[0002] Swell and wind wave are both common wave types in the ocean, the difference between them is that: swell is generally a long-period and long-wavelength wave generated by wind far from the wind source, when strong wind blows the sea water in a large range, swell will also continue to expand outward until the energy is dispersed or encounters an obstacle. Wind wave is a short-period and short-wavelength wave generated by local wind directly acting on the sea surface, which is more local and random. These waves spread over the vast sea surface through wind power, and the effective wave height generated by them affects the safety of navigation, offshore operations and even the shaping of the coastline. The ocean effective wave height can be roughly divided into three states: swell height containing only swell, wind wave height containing only wind wave and mixed wave height containing both swell and wind wave. Since the effective wave height in different states will produce different advantages and disadvantages, it is particularly important to study how to more efficiently and accurately invert the effective wave height in three different states.

[0003] At present, spaceborne GNSS-R technology is widely used in ocean environment research. For ocean effective wave height, the main focus of existing research is to study the wave height of swell and wind wave combination or to establish separate inversion models to invert the effective wave height containing only swell or wind wave. There is no research on a model to invert the effective wave height containing only swell or wind wave, so it is not possible to efficiently and economically completely separate swell height and wind wave height to discuss their respective influences. Considering this aspect, the present application innovatively proposes a method of combining GNSS reflection measurement (GNSS-R) technology and deep learning algorithm to construct a multi-task joint inversion model, hoping to realize the vision of multiple input and output through this method, so as to improve the inversion efficiency and accuracy.

[0004] The construction of a multi-task joint model cannot rule out the occurrence of a tricky nonlinear situation, but deep learning has a very prominent advantage in dealing with nonlinear problems. When constructing the model, the deep learning network can learn the characteristics from simple to complex through a multi-level nonlinear transformation structure, and then select different activation functions (such as ReLU activation function, Sigmoid activation function), loss functions (such as mean square error MSE, Huber loss function HuberLoss) and information propagation methods (such as forward propagation, back propagation) according to the actual situation to further optimize the model. The activation function is the core of the nonlinear capability of the neural network, which can help the neural network to construct a complex nonlinear mapping during modeling. The loss function is generally used to measure the difference between the model output and the true value, and the network weight is adjusted through the back propagation algorithm, which means that the network parameters can be adjusted by optimizing the loss function to obtain better inversion results. During the whole model training process, the information propagation method and gradient update can also be adjusted to further optimize the nonlinear model.

[0005] Considering the above factors and the goal to be achieved by the model, the present application first combines a convolutional neural network (CNN), a dynamic multi-path Mamba model and a Transformer model, on this basis, an adaptive loss layer is added and a regularization loss is introduced, a multi-task joint inversion model capable of realizing multi-input and double-output is constructed, which has almost the same feature capturing ability, generalization ability and flexibility as a general single-output inversion model, and the precision is also comparable, and the efficiency and economic level show the unique advantages of the model. This invention not only first shows the great potential and advantages of a hybrid deep learning model in using multi-modal satellite-borne GNSS-R data to jointly invert ocean swell height and wind wave height, but also promotes the further development of ocean environment research, and provides a thought guidance for the future realization of multi-input and multi-output task model construction.

[0006] In summary, we invented a method for jointly inverting ocean swell height and wind wave height from satellite-borne GNSS-R multi-modal data fusion. SUMMARY

[0007] To solve the above problems, the present application provides a method for jointly inverting ocean swell height and wind wave height from satellite-borne GNSS-R multi-modal data fusion, which is realized by the following technical solutions.

[0008] The method for jointly inverting ocean swell height and wind wave height from satellite-borne GNSS-R multi-modal data fusion comprises the following steps:

[0009] S1, raw data acquisition, acquiring L1 data of Tianmu No. 1 global navigation satellite occultation detector-M, ERA5 reanalysis data set and auxiliary data as raw data;

[0010] S2, pre-processing, quality control and spatio-temporal matching of the original data, and calibration of the L1 data;

[0011] S3, dividing the data set processed in step S2, 60% as training data, 15% as validation data, and 25% as test data;

[0012] S4, constructing a joint inversion model of ocean swell height and wind wave height, and training and validating the model;

[0013] S5, adding an adaptive loss layer, adaptively weighting ocean swell height and wind wave height, introducing regularization to prevent overfitting, and balancing the output of the joint inversion model;

[0014] S6, inputting test data into the validated joint inversion model to further test the fitting effect of the model;

[0015] S7, comparing the output of the joint inversion model with the true value of the swell height and the true value of the wind wave height of the ERA5 data to evaluate the performance of the model.

[0016] Preferably, in step S1,

[0017] The L1 data includes:

[0018] The image-related bistatic radar scattering cross section, effective scattering area, original count DDM and power DDM;

[0019] The DDM-related DDM bistatic radar scattering cross section calculation coefficient, DDM peak signal-to-noise ratio, DDM power calculation coefficient, DDM waveform front slope, DDM specular reflection point normalized bistatic radar scattering cross section, normalized mirror point signal-to-noise ratio, DDM specular reflection point signal-to-noise ratio, DDM sampling intermediate time corresponding to UTC time and DDM waveform data track number;

[0020] The transmitter-related GNSS satellite velocity X component, Y component and Z component corresponding to the signal transmission time of the DDM sampling intermediate time, and the GNSS satellite PRN code used to generate the DDM waveform;

[0021] The receiver-related low-orbit satellite altitude corresponding to the DDM sampling intermediate time, the low-orbit satellite longitude corresponding to the DDM sampling intermediate time, the low-orbit satellite latitude corresponding to the DDM sampling intermediate time, the low-orbit satellite roll angle corresponding to the DDM sampling intermediate time, and the low-orbit satellite velocity X component, Y component and Z component corresponding to the DDM sampling intermediate time;

[0022] height of the specular reflection point, receiver antenna gain of the specular reflection point, azimuth angle of the specular reflection point in the antenna coordinate system, azimuth angle of the specular reflection point in the satellite body coordinate system, azimuth angle of the specular reflection point in the orbit coordinate system, azimuth angle of the specular reflection point in the pattern coordinate system, GNSS signal incidence angle at the specular reflection point, longitude of the specular reflection point, latitude of the specular reflection point, elevation angle of the specular reflection point in the antenna coordinate system, elevation angle of the specular reflection point in the satellite body coordinate system, elevation angle of the specular reflection point in the orbit coordinate system, elevation angle of the specular reflection point in the pattern coordinate system, and X component, Y component and Z component of the velocity of the specular reflection point;

[0023] The ERA5 reanalysis dataset includes wind speed, wind direction, water depth and SWH data;

[0024] The auxiliary data includes ERA5 precipitation data, ERA5 sea surface temperature data, SMAP sea surface salinity data, and OSCAR ocean surface current velocity data along the longitude and latitude directions.

[0025] Preferably, the step S2 comprises the following sub-steps:

[0026] S21, calculating the characteristic observations from the L1 data, performing data quality control using the L1 data product DDM quality identifier, and spatiotemporally matching the L1 data, the ERA5 reanalysis dataset and the auxiliary data;

[0027] S21, mutually calibrating the characteristic observations of the four systems of L1 data:

[0028] First, calibrate the original count DDM to the power DDM, and then calibrate the power DDM to the bistatic radar scattering cross section, which is expressed as follows:

[0029] P = C / G

[0030]

[0031] Where C is the original count DDM, G is the instrument gain, and P is the power DDM, are the bistatic radar scattering cross section, power DDM, and inverse of effective scattering area at delay τ and Doppler f, respectively, P t G t and G r are the equivalent isotropic radiated power of the GNSS transmitter and the receiver antenna gain of the specular reflection point, R t and R r are the distances from the specular reflection point to the transmitter and the receiver, respectively;

[0032] Mutually calibrate the characteristic observations of the four systems of L1 data using a linear model, which is as follows:

[0033] Z calibrated = aZ raw + b

[0034] where Z raw is the uncalibrated raw observation, Z calibrated is the calibrated observation, a is the scale factor to adjust the scale difference between different systems, and b is the bias correction to correct the reference error between different systems.

[0035] Preferably, in the step S4, the main body of the joint inversion model combines the following five input lines:

[0036] The input of the first input line is the bistatic radar cross section;

[0037] The input of the second input line is the GNSS-R variable and auxiliary parameter;

[0038] The input of the third input line is the effective scattering area;

[0039] The input of the fourth input line is the raw count DDM;

[0040] The input of the fifth input line is the power DDM.

[0041] Preferably, the input of the first input line, the fourth input line and the fifth input line are all images, and their processing methods are as follows:

[0042] 1) CNN convolution processing:

[0043] Two convolution layers are used to extract features, and then a global average pooling layer is used to reduce the dimension, and then a fully connected layer is used for further feature mapping and output. This process can be expressed as follows:

[0044] Two convolution layers are used to extract features:

[0045] X1 = Conv2D(X) e R H×W×64

[0046] X2 = Conv2D(X1) e R H×W×128

[0047] where the input tensor is represented by X e R H×W×C , H is the height of the input image, W is the width, and C is the number of convolution kernels; Conv2D() represents a two-dimensional convolution operation; X represents the original input; X1 is the tensor after the first convolution, and the number of convolution kernels is 64; X2 is the tensor after the second convolution, and the number of convolution kernels is 128;

[0048] Enter the global average pooling layer:

[0049] X3 = GlobalAveragePooling2D(X2) e R 128

[0050] Where X3 is a vector of length 128, representing the extracted features, and GlobalAveragePooling2D() represents the global average pooling operation.

[0051] X3 is mapped to 256 neurons through a fully connected layer and outputs X4:

[0052] X4 = Dense(X3) e R 256 ;

[0053] 2) Use the Transformer module to process the input data:

[0054] The Transformer module uses multi-head attention and a feed-forward neural network to process the input data. The attention output is added to the original input through a residual connection, and then layer normalization is applied to the sum to stabilize training and improve convergence. The feed-forward network uses two fully connected layers, the first layer uses a Relu activation function to introduce nonlinearity, and the second layer uses an inputs.shape[-1] operation to project the output to the same dimension as the input. To prevent overfitting, random dropout is used, and the feed-forward output is added to the result of the first layer normalization to apply another residual connection. Finally, the result is normalized again as the final output of the Transformer module. The feed-forward neural network here is implemented by the following formula

[0055] FFN(X) = max(0, XW1 + b1)W2 + b2

[0056] Here W1 represents the first layer weight matrix; W2 represents the second layer weight matrix; b1, b2 represent the bias term; max(0, ·) represents the Relu activation function; X represents the original input;

[0057] The final output can be represented as:

[0058] Transformer(X) = LayerNorm(X + MultiHead(Q, K, V)) + LayerNorm(FFN(X))

[0059] Where MultiHead(Q, K, V) is the multi-head attention mechanism, and LayerNorm() represents the normalization operation.

[0060] Preferably, the input of the third input line is directly processed by a Transformer module, and the four lines of the first input line, the third input line, the fourth input line and the fifth input line are respectively flattened after obtaining the final output of the Transformer module, and then added after passing through a fully connected layer.

[0061] Preferably, the input of the second input line uses a dynamic multi-path state space model Mamba to extract features, uses three different paths to process the input, and each path uses a multi-layer perceptron MLP and a one-dimensional convolution to extract features. The extracted features are combined using multiplication, in addition, a gating mechanism is added, the Softmax gating weight is used to determine the contribution of each path, and the Skip connection combines the final output and the original input to promote the gradient flow. In this module, the input is assumed to be X∈R N×D , N is the batch size, D is the feature dimension, and the final output can be represented by the formula:

[0062]

[0063] where Y final is the weighted sum of each path, is the result of the MLP processing and the one-dimensional convolution Conv1D of each path, and g i represents the gating weight.

[0064] After processing, it passes through a fully connected layer with 128 neurons and uses a ReLU activation function, and is combined with other lines, and the combination process is a data fusion mechanism, which is specifically manifested as multiplication fusion. Before this step, both lines pass through a fully connected layer with 256 neurons and use a ReLU activation function, which is to ensure the consistency of the shape.

[0065] After the five input lines are combined, a bidirectional convolutional long short-term memory network feature reasoning module is further used to effectively extract spatial features in the input data and capture long-term dependencies in the time series data, so as to more comprehensively understand the spatial-temporal features of the data.

[0066] Preferably, step S5 includes the following sub-steps:

[0067] S51, since the final output of the model is two targets, the swell height and the wind wave height, an adaptive loss layer is added to adaptively weight the two targets, which is specifically manifested as follows: first, calculate the loss of each target, then update the weight of each target using adaptive weight, and finally select the target one to pay attention to in the total loss:

[0068]

[0069] where α1 and a 2 is a parameter for adjusting weights, loss1=MSE(y true,1 ,y pred,1 ) and loss2=MSE(y true,2 ,y pred,2 ) are two loss terms;

[0070] S52, regularization is introduced to the weights for constraining a1 and a2 to control the complexity of the model, prevent overfitting, and balance the output of the joint inversion model, then:

[0071]

[0072] The above simultaneously considers the task loss and the regularization loss.

[0073] In the step S7, the performance evaluation indicators of the model are root mean square error RMSE, mean square error MSE, mean absolute error MAE, and Pearson correlation coefficient R 2 For the evaluation indicators of the surge height inversion accuracy, the correlation calculation formula is as follows:

[0074]

[0075] Bias=E[F(x)]-A(x)

[0076]

[0077] Wherein, m is the number of data samples, y i,M is the effective surge height value estimated by the model, y i,T is the surge height value from the reference data set, are the average values of y i,M and y i,T respectively, E[F(x)] is the expectation of the predicted value of the model under a specific input x, A(x) is the true value, A t is the true value at time t, F t is the predicted value at time t, and |·| represents the absolute value operation.

[0078] The wind wave height evaluation indicators are calculated in the same way.

[0079] The beneficial effect of the present invention is that the model of the present invention can directly use this model to output two targets at one time, and to a certain extent, simultaneously ensure the inversion accuracy of the two targets. Compared with various traditional models that can only invert one output target separately, it can greatly save computing efficiency while bringing an effect that is not inferior to the inversion of a single output model, and has economic benefits; it lays the foundation for the multi-input and multi-output design of future models and opens up research ideas in related fields. BRIEF DESCRIPTION OF THE DRAWINGS

[0080] In order to more clearly illustrate the technical solution of the present invention, the following is a brief introduction to the drawings required for use in the description of the specific implementation methods. 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.

[0081] Figure 1 :Flowchart of the joint inversion method of ocean surge height and wind wave height based on multi-modal data fusion of spaceborne GNSS-R;

[0082] Figure 2 : Structure diagram of the joint inversion model of ocean surge height and wind wave height;

[0083] Figure 3 : Scatter density plot of the comparison between the swell height and wind wave height inverted by the dual output of the joint inversion model and the swell height and wind wave height of ERA5;

[0084] Figure 4 : Scatter density plot comparing the swell height and wind wave height retrieved by the joint inversion model with the swell height and wind wave height retrieved by the ERA5;

[0085] Figure 5 : Scatter density plot comparing the swell height and wind wave height retrieved by the single output of DCNN model with the swell height and wind wave height retrieved by ERA5. DETAILED DESCRIPTION

[0086] 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 any creative efforts shall fall within the scope of protection of the present invention.

[0087] Example

[0088] like Figure 1 As shown in FIG, the joint inversion method of ocean surge height and wind wave height by fusion of spaceborne GNSS-R multimodal data includes the following steps:

[0089] S1, original data acquisition, obtaining L1 data of Tianmu No. 1 global navigation satellite occultation detector-M, ERA5 reanalysis data set and auxiliary data as original data;

[0090] S2, pre-processing, quality control and space-time matching of the original data, and calibration of the L1 data;

[0091] S3, dividing the data set processed in step S2, 60% as training data, 15% data as validation data, and 25% as test data;

[0092] S4, constructing a joint inversion model of ocean swell height and wind wave height, and training and validating the model;

[0093] S5, adding an adaptive loss layer, adaptively weighting ocean swell height and wind wave height, introducing regularization to prevent overfitting, and balancing the output of the joint inversion model;

[0094] S6, inputting the test data into the validated joint inversion model to further test the fitting effect of the model;

[0095] S7, comparing the output results of the joint inversion model with the true values of the swell height and the wind wave height of the ERA5 data to evaluate the performance of the model.

[0096] As a further embodiment of the present embodiment, in step S1,

[0097] The L1 data includes:

[0098] The image-related bi-static radar cross section, effective scattering area, original count DDM and power DDM;

[0099] The DDM-related DDM bi-static radar cross section calculation coefficient, DDM peak signal-to-noise ratio, DDM power calculation coefficient, DDM waveform front slope, DDM specular reflection point normalized bi-static radar cross section, normalized mirror point signal-to-noise ratio, DDM specular reflection point signal-to-noise ratio, DDM sampling intermediate time corresponding to UTC time and track number of DDM waveform data;

[0100] The transmitter-related DDM sampling intermediate time corresponding to the GNSS satellite velocity X component, Y component and Z component of the signal transmission time and the GNSS satellite PRN code used to generate the DDM waveform;

[0101] an altitude of the low earth orbit satellite corresponding to the intermediate time of the DDM sampling, a longitude of the low earth orbit satellite corresponding to the intermediate time of the DDM sampling, a latitude of the low earth orbit satellite corresponding to the intermediate time of the DDM sampling, a roll angle of the low earth orbit satellite corresponding to the intermediate time of the DDM sampling, and an X component, a Y component, and a Z component of a velocity of the low earth orbit satellite corresponding to the intermediate time of the DDM sampling;

[0102] an altitude of the specular reflection point, a receiver antenna gain of the specular reflection point, an azimuth angle of the specular reflection point in an antenna coordinate system, an azimuth angle of the specular reflection point in a satellite body coordinate system, an azimuth angle of the specular reflection point in an orbit coordinate system, an azimuth angle of the specular reflection point in a pattern coordinate system, an incident angle of a GNSS signal at the specular reflection point, a longitude of the specular reflection point, a latitude of the specular reflection point, an elevation angle of the specular reflection point in the antenna coordinate system, an elevation angle of the specular reflection point in the satellite body coordinate system, an elevation angle of the specular reflection point in the orbit coordinate system, an elevation angle of the specular reflection point in the pattern coordinate system, and an X component, a Y component, and a Z component of a velocity of the specular reflection point;

[0103] The ERA5 reanalysis dataset includes wind speed, wind direction, water depth, and SWH data;

[0104] The auxiliary data includes ERA5 precipitation data, ERA5 sea surface temperature data, SMAP sea surface salinity data, and OSCAR ocean surface current velocity data along the longitude and latitude directions.

[0105] As a further implementation of the embodiment, step S2 includes the following sub-steps:

[0106] S21, calculating feature observations from the L1 data, performing data quality control using the L1 data product DDM quality identifier, and performing spatio-temporal matching of the L1 data, the ERA5 reanalysis dataset, and the auxiliary data;

[0107] S21, mutually calibrating the feature observations of the four systems of L1 data:

[0108] First, calibrate the original count DDM to the power DDM, and then calibrate the power DDM to the bistatic radar cross section, which is expressed as follows:

[0109] P = C / G

[0110]

[0111] where C is the original count DDM, G is the instrument gain, and P is the power DDM, are the bistatic radar cross section, the power DDM, and the inverse of the effective scattering area at delay τ and Doppler f, respectively, and P t G t and G ris the equivalent isotropically radiated power of the GNSS transmitter and the receiver antenna gain of the specular reflection point, R t and R r are the distances from the specular reflection point to the transmitter and receiver, respectively;

[0112] The four systems of L1 data feature observations are calibrated to each other using a linear model as follows:

[0113] Z calibrated = aZ raw + b

[0114] where Z raw is the uncalibrated raw observation, Z calibrated is the calibrated observation, a is the scale factor to adjust the scale difference between different systems, and b is the bias correction to correct the reference error between different systems.

[0115] As a further implementation of the embodiment, in step S4, the main body of the joint inversion model combines the following five input lines:

[0116] The input of the first input line is the bistatic radar cross section;

[0117] The input of the second input line is the GNSS-R variable and auxiliary parameter;

[0118] The input of the third input line is the effective scattering area;

[0119] The input of the fourth input line is the raw count DDM;

[0120] The input of the fifth input line is the power DDM.

[0121] As a further implementation of the embodiment, as shown in Figure 2 , the inputs of the first input line, the fourth input line, and the fifth input line are all images, and their processing methods are as follows:

[0122] 1) CNN convolution processing:

[0123] Two convolution layers are used to extract features, and then a global average pooling layer is used to reduce the dimension, and then a fully connected layer is used for further feature mapping and output, which can be expressed as follows:

[0124] Two convolution layers are used to extract features:

[0125] X1 = Conv2D(X) e R H×W×64

[0126] X2 = Conv2D(X1) e R H×W×128

[0127] where the input tensor is X e RH×W×C where H is the height and W is the width of the input image, and C is the number of convolution kernels; Conv2D() represents a two-dimensional convolution operation; X represents the original input; X1 is the tensor after the first convolution, and the number of convolution kernels is 64; X2 is the tensor after the second convolution, and the number of convolution kernels is 128;

[0128] Enter the global average pooling layer:

[0129] X3 = GlobalAveragePooling2D(X2) e R 128

[0130] where X3 is a vector with a length of 128, representing the extracted features, and GlobalAveragePooling2D() represents the global average pooling operation;

[0131] Map X3 to 256 neurons through the fully connected layer and output X4:

[0132] X4 = Dense(X3) e R 256 ;

[0133] 2) Use the Transformer module to process the input data:

[0134] The Transformer module uses multi-head attention and a feed-forward neural network to process the input data. The attention output is added to the original input through a residual connection, and then layer normalization is applied to the sum to stabilize training and improve convergence. The feed-forward network uses two fully connected layers. The first layer uses a Relu activation function to introduce nonlinearity, and the second layer uses an inputs.shape[-1] operation to project the output to the same dimension as the input. To prevent overfitting, random dropout is used, and the feed-forward output is added to the result of the first layer normalization to apply another residual connection. Finally, the result is normalized again as the final output of the Transformer module. The feed-forward neural network here is implemented by the following formula

[0135] FFN(X) = max(0, XW1 + b1)W2 + b2

[0136] Here W1 represents the first layer weight matrix; W2 represents the second layer weight matrix; b1, b2 represent the bias term; max(0, ·) represents the Relu activation function; X represents the original input;

[0137] The final output can be represented as:

[0138] Transformer (X) = LayerNorm (X + MultiHead (Q, K, V)) + LayerNorm (FFN (X))

[0139] wherein MultiHead (Q, K, V) is a multi-head attention mechanism, and LayerNorm () represents a normalization operation.

[0140] As a further implementation of the embodiment, the input of the input line three is directly processed by the Transformer module, the input line one, the input line three, the input line four and the input line five are respectively flattened after obtaining the final output of the Transformer module, and then added after passing through the fully connected layer, and merged into one line.

[0141] As a further implementation of the embodiment, the input of the input line two uses a dynamic multi-path state space model Mamba to extract features, uses three different paths to process the input, and each path uses a multi-layer perceptron MLP and a one-dimensional convolution for feature extraction. The extracted features are combined using multiplication, in addition, a gating mechanism is added, the Softmax gating weight is used to determine the contribution of each path, and the Skip connection combines the final output and the original input to promote gradient flow. In this module, the input is assumed to be X N×D , N is the batch size, D is the feature dimension, and the final output can be represented by the formula:

[0142]

[0143] wherein Y final is the weighted sum of each path, is the result of MLP processing and one-dimensional convolution Conv1D through each path, and g i represents the gating weight.

[0144] After processing, it passes through a fully connected layer with 128 neurons and uses a ReLU activation function, and is combined with other lines, and the combination process is a data fusion mechanism, which is specifically manifested as multiplication fusion. Before this step, both lines pass through a fully connected layer with 256 neurons and use a ReLU activation function, which is to ensure the consistency of the shape.

[0145] After the five input lines are combined, a bidirectional convolutional long short-term memory network feature reasoning module is further used to effectively extract spatial features in the input data and capture long-term dependencies in the time series data, so as to more comprehensively understand the space-time features of the data.

[0146] As a further implementation of the embodiment, step S5 comprises the following sub-steps:

[0147] S51, since the model finally outputs two targets, respectively, the swell height and the wind wave height, an adaptive loss layer is added to adaptively weight the two targets, which is specifically manifested as follows: first, the loss of each target is calculated, then the adaptive weight is used to update the weight of each target, and finally the total loss selects the attention of the first target:

[0148]

[0149] Wherein, α 1 and α 2 are parameters for adjusting weights, loss1=MSE(y true,1 ,y pred,1 ) and loss2=MSE(y true,2 ,y pred,2 ) are two loss terms.

[0150] S52, the weight is introduced into the regularization for constraining α1 and α2, so as to control the complexity of the model, prevent overfitting, and balance the output of the joint inversion model, then:

[0151]

[0152] The above considers two aspects of task loss and regularization loss.

[0153] In step S7, the performance evaluation index of the model is root mean square error RMSE, mean square error MSE, mean absolute error MAE and Pearson correlation coefficient R 2 For the swell height inversion accuracy evaluation index, the correlation calculation formula is as follows:

[0154]

[0155] Bias=E[F(x)]-A(x)

[0156]

[0157] Wherein, m is the number of data samples, y i,M is the effective swell height value estimated by the model, y i,T is the swell height value from the reference data set, are the average values of y i,M and y i,T , E[F(x)] is the expectation of the predicted value of the model under a specific input x, A(x) is the true value, A t is the true value at time t, F t is the predicted value at time t, and |·| represents the absolute value operation.

[0158] The wind wave height evaluation index is calculated in the same way as above.

[0159] One specific embodiment of the present application is as follows:

[0160] The L1 data of the Tianmu No. 1 global navigation satellite occultation detector-M is provided by Space Tianmu (Chongqing) Satellite Technology Co., Ltd. and Space Technology Co., Ltd. (Beijing) Space Information Application Co., Ltd.

[0161] The data used covers the period from January to July 2024.

[0162] The basic configuration of the experimental platform for constructing the ocean swell height and wind wave height joint inversion model of Tianmu No. 1 satellite-borne GNSS-R multi-modal data fusion is shown in Table 1:

[0163]

[0164] Two experiments were conducted using the ocean swell height and wind wave height joint inversion model, and two experimental results were given. One is the joint inversion model double output inversion swell height and wind wave height result, and the other is the joint inversion model single output inversion swell height and single output inversion wind wave height result. In addition to comparing the experimental results of the two with each other, a deep learning method (DCNN) is used for comparison and verification. Four performance evaluation indexes are calculated for comparison with ERA5 data, and the experimental results are shown in Table 2:

[0165] Table 2 Comparison of precision statistics of three experimental results with ERA5 data

[0166]

[0167] From the above table, the swell height mean square error RMSE of the joint inversion model double output is 0.586, which is reduced by 1.72% and 9.89% compared with the joint inversion model single output and the DCNN model single output, respectively. The wind wave height mean square error RMSE of the joint inversion model double output is 0.330, which is reduced by 3.28% and 0.94% compared with the joint inversion model single output and the DCNN model single output, respectively. In addition, the comparison of scatter point density graphs of the joint inversion model double output inversion swell height and wind wave height with ERA5 swell height and wind wave height, the comparison of scatter point density graphs of the joint inversion model single output inversion swell height and single output inversion wind wave height with ERA5 swell height and wind wave height, and the comparison of scatter point density graphs of the single output inversion swell height and single output inversion wind wave height using the DCNN model with ERA5 swell height and wind wave height are obtained, as shown in Figures 3-5From the experimental results, the joint inversion model double output of the application has a slight advantage in inversion accuracy compared with the joint inversion model single output and the DCNN model single output, but the efficiency and economy highlight the unique advantages of the model.

[0168] Compared with many GNSS-R sea surface significant wave height inversion methods combined with deep learning, the advantage of the model of the application is that, unlike various models in the past that can only individually invert one output target, the model can directly output two targets at once, and to some extent, ensures the inversion accuracy of the two targets. This method greatly saves the calculation efficiency while bringing the same effect as the single output model inversion, and when applied to the economic level, it also saves the cost accordingly.

[0169] In summary, the implementation of the application is not limited to the present embodiment, but also lays the foundation for the design of multiple input and multiple output models in the future, and opens up research ideas in related fields. This also shows that the practical significance of the spaceborne GNSS-R technology to invert the significant wave height is continuously strengthening.

[0170] The preferred embodiments of the application disclosed above are only used to help explain the application. The preferred embodiments do not describe all the details and limit the application to the specific embodiments. Obviously, many modifications and changes can be made according to the content of the specification. The specification selects and describes these embodiments in order to better explain the principles and practical applications of the application, so that those skilled in the art can well understand and use the application. The application is limited by the claims and their entire scope and equivalents.

Claims

1. A method for ocean swell height and wind wave height joint inversion of space-borne GNSS-R multi-modal data fusion, characterized in that, The method comprises the following steps: S1, original data acquisition, obtaining L1 data of Tianmu No. 1 global navigation satellite occultation detector-M, ERA5 reanalysis data set and auxiliary data as original data; S2, preprocessing, quality control and space-time matching of the original data, and calibration of the L1 data; S3, dividing the data set processed in step S2, 60% as training data, 15% as verification data, and 25% as test data; S4, constructing an ocean swell height and wind wave height joint inversion model, and training and verifying the model; The main body of the joint inversion model combines the following five input lines: The input of the input line one is the bistatic radar scattering cross section; The input of the input line two is the GNSS-R variable and auxiliary parameter; The input of the input line three is the effective scattering area; The input of the input line four is the original count DDM; The input of the input line five is the power DDM; The inputs of the input line one, the input line four and the input line five are all images, and their processing methods are as follows: 1) CNN convolution processing: Two convolution layers are used to extract features, then a global average pooling layer is used to reduce the dimension, and then a fully connected layer is used for further feature mapping and output, which is expressed by the following formula: Two convolution layers are used to extract features: X1= Conv2D(X) e R H×W×64 X2 = Conv2D(X1) e R H×W×128 wherein the input tensor is denoted by X ∈ R H×W×C H is the height of the input image, W is the width, and C is the number of convolution kernels; Conv2D() represents a two-dimensional convolution operation; X represents the original input; X1 is a tensor after the first convolution, and the number of convolution kernels is 64; X2 is a tensor after the second convolution, and the number of convolution kernels is 128; Enter the global average pooling layer: X3 = GlobalAveragePooling2D(X2) e R 128 Wherein, X3 is a vector with a length of 128, representing the extracted features, and GlobalAveragePooling2D() represents the global average pooling operation; X4 is mapped to 256 neurons through the fully connected layer and output: X4 = Dense(X3) e R 256 ; 2) using the Transformer module to process the input data: The Transformer module uses multi-head attention and a feedforward neural network to process the input data, the attention output is added to the original input through a residual connection, then the layer normalization is applied to the sum to stabilize the training and improve the convergence effect, the feedforward network uses two fully connected layers, the first layer uses the Relu activation function to introduce nonlinearity, the second layer makes inputs.shape[-1] operation, and the output is projected to the same dimension as the input, in order to prevent overfitting, random deactivation is used, and the feedforward output is added to the result of the first layer normalization to apply another residual connection, and the result is normalized again as the final output of the Transformer module, the feedforward neural network here is realized by the following formula FFN(X)=max(0,XW1+b1)W2+b2 Here W1 represents the first layer weight matrix; W2 represents the second layer weight matrix; b1, b2 represent the bias term; max(0, ·) represents the Relu activation function; X represents the original input; The final output can be represented as: Transformer(X)=LayerNorm(X+MultiHead(Q,K,V))+LayerNorm(FFN(X)) Wherein, MultiHead(Q,K,V) is a multi-head attention mechanism, and LayerNorm() represents the normalization operation; The input of the second input line uses a dynamic multi-path state space model Mamba to extract features, uses three different paths to process the input, and each path uses a multi-layer perceptron MLP and a one-dimensional convolution for feature extraction. The extracted features are combined using multiplication, and in addition, a gating mechanism is added to determine the contribution of each path using Softmax gate weights. The Skip connection combines the final output and the original input, facilitating gradient flow. In this module, the input is assumed to be X∈R N ×D , N is the batch size, and D is the feature dimension. The final output can be represented by the formula: where Y final is the weighted sum of each path, is the result of the MLP processing and one-dimensional convolution Conv1D through each path, g i denotes the gating weight; After processing, the two lines are merged through a fully connected layer with 128 neurons and using ReLU activation function, and the merging process is the data fusion mechanism, which is specifically manifested as multiplication fusion. Before this step, both lines have passed through a fully connected layer with 256 neurons and using ReLU activation function, which is to ensure the consistency of the shape; After the five input lines are merged, a bidirectional convolutional long short-term memory network feature reasoning module is further used to effectively extract the spatial features in the input data and capture the long-term dependence in the time series data, so as to more comprehensively understand the spatial-time features of the data; S5, an adaptive loss layer is added to adaptively weight the two targets of ocean swell height and wind wave height, and regularization is introduced to the weight to prevent overfitting and balance the output of the joint inversion model; S6, the test data is input into the verified joint inversion model to further test the fitting effect of the model; S7, the output results of the joint inversion model are compared with the true values of the swell height and the wind wave height of the ERA5 data to evaluate the performance of the model.

2. The space-borne GNSS-R multi-modal data fusion combined ocean swell and wind wave height retrieval method according to claim 1, characterized in that, In the step S1, The L1 data includes: The image-related bistatic radar scattering cross section, effective scattering area, original count DDM and power DDM; The DDM-related DDM bistatic radar scattering cross section calculation coefficient, DDM peak signal-to-noise ratio, DDM power calculation coefficient, DDM waveform front slope, DDM specular reflection point normalized bistatic radar scattering cross section, normalized mirror point signal-to-noise ratio, DDM specular reflection point signal-to-noise ratio, DDM sampling intermediate time corresponding UTC time and DDM waveform data track number; The transmitter-related GNSS satellite velocity X component, Y component and Z component corresponding to the signal transmission time of the DDM sampling intermediate time and the GNSS satellite PRN code used to generate the DDM waveform; The receiver-related low-orbit satellite height corresponding to the DDM sampling intermediate time, low-orbit satellite longitude corresponding to the DDM sampling intermediate time, low-orbit satellite latitude corresponding to the DDM sampling intermediate time, low-orbit satellite roll angle corresponding to the DDM sampling intermediate time, and low-orbit satellite velocity X component, Y component and Z component corresponding to the DDM sampling intermediate time; The specular reflection point height related to the geometric relationship, specular reflection point receiver antenna gain, specular reflection point azimuth in the antenna coordinate system, specular reflection point azimuth in the satellite body coordinate system, specular reflection point azimuth in the orbit coordinate system, specular reflection point azimuth in the directional diagram coordinate system, specular reflection point GNSS signal incident angle, specular reflection point longitude, specular reflection point latitude, specular reflection point elevation angle in the antenna coordinate system, specular reflection point elevation angle in the satellite body coordinate system, specular reflection point elevation angle in the orbit coordinate system, specular reflection point elevation angle in the directional diagram coordinate system, and specular reflection point velocity X component, Y component and Z component; The ERA5 reanalysis dataset includes wind speed, wind direction, water depth and SWH data; The auxiliary data includes ERA5 precipitation data, ERA5 sea surface temperature data, SMAP sea surface salinity data, and OSCAR ocean surface flow velocity data along the longitude and latitude directions.

3. The method for joint inversion of ocean surge height and wind wave height by fusion of spaceborne GNSS-R multimodal data according to claim 2 is characterized in that: The step S2 comprises the following sub-steps: S21, calculating feature observations from L1 data, using L1 data product DDM quality identifier for data quality control, spatio-temporally matching L1 data, ERA5 reanalysis dataset and auxiliary data; S21, mutually calibrating feature observations of four systems of L1 data: Firstly, calibrate the original count DDM to the power DDM, and then calibrate the power DDM to the bistatic radar scattering cross section, which is expressed as follows: where C is the raw count DDM, G is the instrument gain, and P is the power DDM, respectively, the bistatic radar cross section, power DDM, and inverse of the effective scattering area at delay τ and Doppler f, P t G t and G r are the GNSS transmitter equivalent isotropically radiated power and the receiver antenna gain of the specular reflection point, R t and R r are the distances from the specular reflection point to the transmitter and receiver, respectively; Mutually calibrate feature observations of four systems of L1 data using a linear model, which is as follows: Z calibrated = αZ raw + β where Z raw is the uncalibrated raw observation, Z calibrated is the calibrated observation, a is a scale factor to adjust the scale difference between different systems, and β is a bias correction to correct the reference error between different systems.

4. The space-borne GNSS-R multi-modal data fusion method for ocean swell and wind wave height joint inversion according to claim 3, characterized in that, The input line three is directly processed by using the Transformer module, and the input line one, the input line three, the input line four and the input line five are respectively flattened after obtaining the final output of the Transformer module, and then added after passing through the full connection layer, and finally combined into one line.

5. The space-borne GNSS-R multi-modal data fusion-based ocean swells and wind waves height joint inversion method according to claim 4, characterized in that, The step S5 comprises the following sub-steps: S51, since the final output of the model is two targets, namely the swell height and the wind wave height, an adaptive loss layer is added to adaptively weight the two targets, which is specifically as follows: firstly, calculate the loss of each target, then update the weight of each target by using the adaptive weight, and finally, the total loss selects the target one to be focused on: wherein, α 1 and α 2 are parameters of the adjustment weight, loss1 = MSE(y true,1 , y pred,1 ) and loss2 = MSE(y true,2 , y pred,2 ) are two loss terms; S52, regularization is introduced on the weights To constrain α1 and α2, so as to control the complexity of the model, prevent overfitting, and balance the output of the joint inversion model, then: The above simultaneously considers the task loss and the regularization loss.

6. The space-borne GNSS-R multi-modal data fusion-based ocean swells and wind waves height joint inversion method according to claim 5, characterized in that, In the step S7, the performance evaluation indexes of the model are root mean square error RMSE, mean square error MSE, mean absolute error MAE and Pearson correlation coefficient R 2 For the evaluation index of the surge height inversion accuracy, the correlation calculation formula is as follows: where m is the number of data samples, y i,M is the model-estimated significant wave height value, y i,T is the significant wave height value from the reference dataset, are the mean values of y i,M and y i,T respectively, E[F(x)] is the expectation of the model’s prediction value at a particular input x, A(x) is the true value, A t is the true value at time t, F t is the predicted value at time t, and |·| denotes the absolute value operation. The wind wave height evaluation index is calculated in the same way as above.

Images