Ocean surge height and wind wave height joint inversion method based on satellite-borne GNSS-R multi-modal data fusion

Abstract

The invention discloses a satellite-borne GNSS-R multi-modal data fusion ocean surge height and wind wave height joint inversion method, and belongs to the field of GNSS reflection measurement and deep learning cross research, and the method 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 an ocean surge height and wind wave height joint inversion model; s5, adaptive weighting is carried out on the two targets of the sea surge height and the wind wave height; s6, inputting the test data into the joint inversion model for testing; the method has the advantages that the model can be directly used for outputting two targets at a time, the inversion accuracy of the two targets is guaranteed to a certain degree at the same time, the calculation efficiency is greatly saved, meanwhile, the effect which is not inferior to that of single-output model inversion can be achieved, and the method is suitable for large-scale popularization and application. And economic benefits are achieved.

Description

Technical Field

[0001] The present invention relates to the interdisciplinary research field of GNSS reflection measurement and deep learning, and in particular to a joint inversion method for ocean surge height and wind wave height by fusing satellite-borne GNSS-R multimodal data. Background Art

[0002] Swells and wind waves are both common types of waves in the ocean. The difference is that swells are generally waves with longer periods and wavelengths generated by wind far away from the wind source. When strong winds continue to blow the seawater over a large area, swells will continue to expand outward until the energy is dispersed or encounter obstacles. Wind waves are waves with shorter periods and wavelengths generated by local wind directly acting on the ocean surface, which are more local and random. These waves are propagated on the vast sea surface by wind power, and the effective wave heights generated affect navigation safety, marine operations and even the shaping of the coastline. The effective wave height of the ocean can be roughly divided into three states: swell height containing only swells, wind wave height containing only wind waves, and mixed wave height containing both swells and wind waves. Since the effective wave heights in different states will have different advantages and disadvantages, it is particularly important to study how to more efficiently and accurately invert the effective wave heights in three different states.

[0003] At present, satellite-borne GNSS-R technology is widely used in marine environmental research. For the effective wave height of the ocean, the main focus of existing research is to study the wave height of the combination of swell waves and wind waves or to establish inversion models to separately invert the effective wave height containing only swell waves or only wind waves. There is no model to invert the effective wave height containing only swell waves or only wind waves. Therefore, it is not possible to efficiently and economically separate the swell height and wind wave height to discuss their respective impacts. Taking this aspect into consideration, the present invention innovatively proposes a method of combining GNSS reflection measurement (GNSS-R) technology with a deep learning algorithm to construct a multi-task joint inversion model, hoping to achieve the vision of multi-input and multi-output through this method, thereby improving the inversion efficiency and accuracy.

[0004] Building a multi-task joint model cannot rule out the possibility of thorny nonlinear situations, but deep learning happens to have a very prominent advantage in dealing with nonlinear problems. When building a model, the deep learning network can learn features 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, backpropagation) according to actual conditions to further optimize the model. The activation function is the core of the nonlinear ability of the neural network, and can help the neural network build complex nonlinear mappings when modeling. The loss function is generally used to measure the difference between the model output and the true value, and adjust the network weights through the backpropagation algorithm, which means that the overall parameters of the network can be adjusted by optimizing the loss function, thereby obtaining better inversion results. During the entire model training process, the information propagation method and gradient update can also be adjusted to further optimize the nonlinear model.

[0005] Taking the above factors into consideration and the goals that the model is expected to achieve, the present invention combines the convolutional neural network (CNN), the dynamic multi-path Mamba model and the Transformer model for the first time, and on this basis, adds an adaptive loss layer and introduces a regularization loss in a targeted manner to construct a multi-task joint inversion model that can achieve multi-input and dual-output. While achieving almost the same feature capture ability, generalization ability and flexibility as the general single-output inversion model, the accuracy is also comparable, and the efficiency and economy show the unique advantages of the model. This invention not only demonstrates for the first time the great potential and advantages of the hybrid deep learning model in using multi-modal satellite-borne GNSS-R data to jointly invert ocean surge height and wind wave height, but also promotes the further development of marine environment research, and at the same time provides ideological guidance for the future construction of multi-input and multi-output task models.

[0006] In summary, we invented a joint inversion method for ocean surge height and wind wave height by fusing spaceborne GNSS-R multimodal data. Summary of the invention

[0007] To solve the above problems, the present invention provides a method for joint inversion of ocean surge height and wind wave height by fusion of satellite-borne GNSS-R multi-modal data, which is achieved through the following technical solutions.

[0008] The joint inversion method of ocean surge height and wind wave height based on spaceborne GNSS-R multi-modal data fusion includes the following steps:

[0009] S1, raw data acquisition, obtains the L1 data of the Tianmu-1 global navigation satellite occultation sounder-M type, the ERA5 reanalysis data set and auxiliary data as raw data;

[0010] S2, preprocessing, quality control and spatiotemporal matching of raw data, and calibration of L1 data;

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

[0012] S4, construct a joint inversion model of ocean surge height and wind wave height, and train and verify the model;

[0013] S5, adds an adaptive loss layer to adaptively weight the two targets of ocean surge height and wind wave height, introduces regularization on the weights to prevent overfitting and balance the output of the joint inversion model;

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

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

[0016] Preferably, in step S1,

[0017] L1 data includes:

[0018] Bistatic radar cross section, effective scattering area, raw count DDM and power DDM associated with the image;

[0019] DDM related DDM bistatic radar cross section calculation coefficient, DDM peak signal-to-noise ratio, DDM power calculation coefficient, DDM waveform leading edge slope, DDM mirror reflection point normalized bistatic radar cross section, normalized mirror point signal-to-noise ratio, DDM mirror reflection point signal-to-noise ratio, UTC time corresponding to the middle moment of DDM acquisition and trajectory number of DDM waveform data;

[0020] The GNSS satellite velocity X component, Y component and Z component at the signal transmission time corresponding to the DDM acquisition intermediate time associated with the transmitter and the GNSS satellite PRN code used to generate the DDM waveform;

[0021] The altitude of the low-orbit satellite corresponding to the middle moment of DDM acquisition associated with the receiver, the longitude and latitude of the low-orbit satellite corresponding to the middle moment of DDM acquisition, the roll angle of the low-orbit satellite corresponding to the middle moment of DDM acquisition, and the X component, Y component, and Z component of the velocity of the low-orbit satellite corresponding to the middle moment of DDM acquisition;

[0022] The height of the mirror reflection point related to the geometric relationship, the gain of the receiver antenna at the mirror reflection point, the azimuth of the mirror reflection point in the antenna coordinate system, the azimuth of the mirror reflection point in the satellite body coordinate system, the azimuth of the mirror reflection point in the orbit coordinate system, the azimuth of the mirror reflection point in the pattern coordinate system, the incident angle of the GNSS signal at the mirror reflection point, the longitude and latitude of the mirror reflection point, the altitude angle of the mirror reflection point in the antenna coordinate system, the altitude angle of the mirror reflection point in the satellite body coordinate system, the altitude angle of the mirror reflection point in the orbit coordinate system, the altitude angle of the mirror reflection point in the pattern coordinate system, and the X component, Y component and Z component of the velocity of the mirror reflection point;

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

[0024] Auxiliary data include ERA5 precipitation data, ERA5 sea surface temperature data, SMAP sea surface salinity data, and OSCAR ocean surface current velocity data along longitude and latitude.

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

[0026] S21, calculate characteristic observations from L1 data, perform data quality control using the L1 data product DDM quality identifier, and perform spatiotemporal matching of L1 data, ERA5 reanalysis datasets, and auxiliary data;

[0027] S21, mutual calibration of characteristic observations of the four systems of L1 data:

[0028] First, the raw count DDM is calibrated to the power DDM, and then the power DDM is calibrated to the bistatic radar cross section. The formula is as follows:

[0029] P=C / G

[0030]

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

[0032] The characteristic observations of the four systems of L1 data are calibrated with each other using a linear model. The linear model is as follows:

[0033] Z calibrated =αZ raw +β

[0034] Among them, Z raw is the uncalibrated raw observation, Z calibrated is the calibrated observation value, α is the scale factor used to adjust the scale difference between different systems, and β is the bias correction used to correct the benchmark error between different systems.

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

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

[0037] The input of input line 2 is GNSS-R variables and auxiliary parameters;

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

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

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

[0041] Preferably, the inputs of input line 1, input line 4 and input line 5 are all images, and their processing methods are as follows:

[0042] 1) CNN convolution processing:

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

[0044] Use two convolutional layers to extract features:

[0045] X 1 =Conv2D(X)∈R H×W×64

[0046] X 2 =Conv2D(X 1 )∈R H×W×128

[0047] Among them, the input tensor is denoted by X∈R H×W×C Indicates that 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; X 1 is the tensor after the first convolution, with 64 convolution kernels; X 2 It is the tensor after the second convolution, and the number of convolution kernels is 128;

[0048] Enter the global average pooling layer:

[0049] X 3 =GlobalAveragePooling2D(X 2 )∈R 128

[0050] Among them, X 3 is a vector of length 128, representing the extracted features, and GlobalAveragePooling2D() represents the global average pooling operation;

[0051] Through the fully connected layer, X 3 Mapped to 256 neurons and output X 4 :

[0052] X 4 =Dense(X 3 )∈R 256 ;

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

[0054] The Transformer module uses multi-head attention and feedforward neural networks to process 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 achieve stable training and improve convergence. The feedforward network uses two fully connected layers. The first layer uses the Relu activation function to introduce nonlinearity. The second layer uses inputs.shape[-1] to project the output to the same dimension as the input. In order to prevent overfitting, random inactivation is used, and the feedforward output is added to the normalized result of the first layer to apply another residual connection. Finally, the result is normalized again as the final output of the Transformer module. The feedforward neural network here is implemented by the following formula

[0055] FFN(X)=max(0,XW 1 +b 1 )W 2 +b 2

[0056] Here W 1 represents the first layer weight matrix; W 2 represents the second layer weight matrix; b 1 , b 2 represents the bias term; max(0, ) represents the Relu activation function; X represents the original input;

[0057] The final output can be expressed as:

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

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

[0060] Preferably, the input of input line three is directly processed by the Transformer module, and the four lines of input line one, input line three, input line four and input line five are flattened after obtaining the final output of the Transformer module, and are added after passing through the fully connected layer to form one line.

[0061] Preferably, the input of input line 2 uses the 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, and uses multiplication combination for the extracted features. In addition, a gating mechanism is added, and the Softmax gating weight is used to determine the contribution of each path, and the Skip connection combines the final output with the original input to promote gradient flow. In this module, assume that the input is X∈R N×D , N is the batch size, D is the feature dimension, and the final output can be expressed by the formula:

[0062]

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

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

[0065] After the five input lines are merged, a bidirectional convolutional long short-term memory network feature inference module is further used to effectively extract the spatial features in the input data while capturing the long-term dependencies in the time series data, thereby more comprehensively understanding the spatial-temporal characteristics of the data.

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

[0067] S51, since the final output of the model is two targets, namely, swell height and wind wave height, an adaptive loss layer is added to adaptively weight the two targets. Specifically, the losses of the two targets are first calculated, and then the adaptive weights are used to update their respective weights. Finally, the total loss is selected to strengthen the focus on target one:

[0068]

[0069] Among them, α 1 and α 2 To adjust the weight parameters, loss 1 =MSE(y true,1 ,y pred,1 ) and loss 2 =MSE(y true,2 ,y pred,2 ) are two loss terms;

[0070] S52, regularization is introduced on weights To constrain α 1 and α 2 , in order to control the complexity of the model, prevent overfitting, and balance the output of the joint inversion model, then:

[0071]

[0072] The above considers both task loss and regularization loss.

[0073] In 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 index of high surge inversion accuracy, the relevant calculation formula is as follows:

[0074]

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

[0076]

[0077] Where m is the number of data samples, y i,M is the effective surge height estimated by the model, y i,T is the swell height value from the reference dataset, y i,M and i,T The average value of E[F(x)] is the expected value of the model for a specific input x, A(x) is the true value, and A t is the true value at time t, F tis the predicted value at time t, and |·| represents an absolute value operation.

[0078] The calculation method of wind and wave height assessment index is the same as above.

[0079] The beneficial effect of the present invention is that the model of the present invention can directly use the present 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 the calculation efficiency while bringing the effect that is not inferior to the inversion of the single output model, and has economic benefits; it lays the foundation for the multi-input and multi-output design of the model in the future, 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 drawings required for use in the description of the specific implementation methods will be briefly introduced below. 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 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 diagram 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 diagram of the comparison between the single output inversion swell height and single output inversion wind wave height of the joint inversion model and the swell height and wind wave height of ERA5;

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

[0086] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. 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 creative work are within the scope of protection of the present invention.

[0087] Example

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

[0089] S1, raw data acquisition, obtains the L1 data of the Tianmu-1 global navigation satellite occultation sounder-M type, the ERA5 reanalysis data set and auxiliary data as raw data;

[0090] S2, preprocessing, quality control and spatiotemporal matching of raw data, and calibration of L1 data;

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

[0092] S4, construct a joint inversion model of ocean surge height and wind wave height, and train and verify the model;

[0093] S5, adds an adaptive loss layer to adaptively weight the two targets of ocean surge height and wind wave height, introduces regularization on the weights to prevent overfitting and balance the output of the joint inversion model;

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

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

[0096] As a further implementation of this embodiment, in step S1,

[0097] L1 data includes:

[0098] Bistatic radar cross section, effective scattering area, raw count DDM and power DDM associated with the image;

[0099] DDM related DDM bistatic radar cross section calculation coefficient, DDM peak signal-to-noise ratio, DDM power calculation coefficient, DDM waveform leading edge slope, DDM mirror reflection point normalized bistatic radar cross section, normalized mirror point signal-to-noise ratio, DDM mirror reflection point signal-to-noise ratio, UTC time corresponding to the middle moment of DDM acquisition and trajectory number of DDM waveform data;

[0100] The GNSS satellite velocity X component, Y component and Z component at the signal transmission time corresponding to the DDM acquisition intermediate time associated with the transmitter and the GNSS satellite PRN code used to generate the DDM waveform;

[0101] The altitude of the low-orbit satellite corresponding to the middle moment of DDM acquisition associated with the receiver, the longitude and latitude of the low-orbit satellite corresponding to the middle moment of DDM acquisition, the roll angle of the low-orbit satellite corresponding to the middle moment of DDM acquisition, and the X component, Y component, and Z component of the velocity of the low-orbit satellite corresponding to the middle moment of DDM acquisition;

[0102] The height of the mirror reflection point related to the geometric relationship, the gain of the receiver antenna at the mirror reflection point, the azimuth of the mirror reflection point in the antenna coordinate system, the azimuth of the mirror reflection point in the satellite body coordinate system, the azimuth of the mirror reflection point in the orbit coordinate system, the azimuth of the mirror reflection point in the pattern coordinate system, the incident angle of the GNSS signal at the mirror reflection point, the longitude and latitude of the mirror reflection point, the altitude angle of the mirror reflection point in the antenna coordinate system, the altitude angle of the mirror reflection point in the satellite body coordinate system, the altitude angle of the mirror reflection point in the orbit coordinate system, the altitude angle of the mirror reflection point in the pattern coordinate system, and the X component, Y component and Z component of the velocity of the mirror reflection point;

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

[0104] Auxiliary data include ERA5 precipitation data, ERA5 sea surface temperature data, SMAP sea surface salinity data, and OSCAR ocean surface current velocity data along longitude and latitude.

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

[0106] S21, calculate characteristic observations from L1 data, perform data quality control using the L1 data product DDM quality identifier, and perform spatiotemporal matching of L1 data, ERA5 reanalysis datasets, and auxiliary data;

[0107] S21, mutual calibration of characteristic observations of the four systems of L1 data:

[0108] First, the raw count DDM is calibrated to the power DDM, and then the power DDM is calibrated to the bistatic radar cross section. The formula is as follows:

[0109] P=C / G

[0110]

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

[0112] The characteristic observations of the four systems of L1 data are calibrated with each other using a linear model. The linear model is as follows:

[0113] Z calibrated =αZ raw +β

[0114] Among them, Z raw is the uncalibrated raw observation, Z calibrated is the calibrated observation value, α is the scale factor used to adjust the scale difference between different systems, and β is the bias correction used to correct the benchmark error between different systems.

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

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

[0117] The input of input line 2 is GNSS-R variables and auxiliary parameters;

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

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

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

[0121] As a further implementation method of this embodiment, Figure 2 As shown, the inputs of input line 1, input line 4 and input line 5 are all images, and their processing methods are as follows:

[0122] 1) CNN convolution processing:

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

[0124] Use two convolutional layers to extract features:

[0125] X 1 =Conv2D(X)∈R H×W×64

[0126] X 2 =Conv2D(X 1 )∈R H×W×128

[0127] Among them, the input tensor is denoted by X∈R H×W×C Indicates that 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; X 1 is the tensor after the first convolution, with 64 convolution kernels; X 2 It is the tensor after the second convolution, and the number of convolution kernels is 128;

[0128] Enter the global average pooling layer:

[0129] X 3 =GlobalAveragePooling2D(X 2 )∈R 128

[0130] Among them, X 3 is a vector of length 128, representing the extracted features, and GlobalAveragePooling2D() represents the global average pooling operation;

[0131] Through the fully connected layer, X 3 Mapped to 256 neurons and output X 4 :

[0132] X 4 =Dense(X 3 )∈R 256 ;

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

[0134] The Transformer module uses multi-head attention and feedforward neural networks to process 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 achieve stable training and improve convergence. The feedforward network uses two fully connected layers. The first layer uses the Relu activation function to introduce nonlinearity. The second layer uses inputs.shape[-1] to project the output to the same dimension as the input. In order to prevent overfitting, random inactivation is used, and the feedforward output is added to the normalized result of the first layer to apply another residual connection. Finally, the result is normalized again as the final output of the Transformer module. The feedforward neural network here is implemented by the following formula

[0135] FFN(X)=max(0,XW 1 +b 1 )W 2 +b 2

[0136] Here W 1 represents the first layer weight matrix; W 2 represents the second layer weight matrix; b 1 , b 2 represents the bias term; max(0, ) represents the Relu activation function; X represents the original input;

[0137] The final output can be expressed as:

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

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

[0140] As a further implementation method of this embodiment, the input of input line three is directly processed by the Transformer module. After obtaining the final output of the Transformer module, the four lines of input line one, input line three, input line four and input line five are flattened respectively, passed through the fully connected layer respectively, and then added together to form one line.

[0141] As a further implementation method of this embodiment, the input of input line 2 uses the 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, and uses multiplication combination for the extracted features. In addition, a gating mechanism is added, and the contribution of each path is determined by the Softmax gating weight, and the Skip connection combines the final output with the original input to promote gradient flow. In this module, assume that the input is X∈R N×D , N is the batch size, D is the feature dimension, and the final output can be expressed by the formula:

[0142]

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

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

[0145] After the five input lines are merged, a bidirectional convolutional long short-term memory network feature inference module is further used to effectively extract the spatial features in the input data while capturing the long-term dependencies in the time series data, thereby more comprehensively understanding the spatial-temporal characteristics of the data.

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

[0147] S51, since the final output of the model is two targets, namely, swell height and wind wave height, an adaptive loss layer is added to adaptively weight the two targets. Specifically, the losses of the two targets are first calculated, and then the adaptive weights are used to update their respective weights. Finally, the total loss is selected to strengthen the focus on target one:

[0148]

[0149] Among them, α 1 and α 2 To adjust the weight parameters, loss 1 =MSE(y true,1 ,y pred,1 ) and loss 2 =MSE(y true,2 ,y pred,2 ) are two loss terms;

[0150] S52, regularization is introduced on weights To constrain α 1 and α 2 , in order to control the complexity of the model, prevent overfitting, and balance the output of the joint inversion model, then:

[0151]

[0152] The above considers both task loss and regularization loss.

[0153] In 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 index of high surge inversion accuracy, the relevant calculation formula is as follows:

[0154]

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

[0156]

[0157] Where m is the number of data samples, y i,M is the effective surge height estimated by the model, y i,T is the swell height value from the reference dataset, y i,M and i,T The average value of E[F(x)] is the expected value of the model for a specific input x, A(x) is the true value, and A t is the true value at time t, F t is the predicted value at time t, and |·| represents an absolute value operation.

[0158] The calculation method of wind and wave height assessment index is the same as above.

[0159] A specific embodiment of the present invention is as follows:

[0160] The L1 data of the Tianmu-1 global navigation satellite occultation detector-M type is provided by Aerospace Tianmu (Chongqing) Satellite Technology Co., Ltd. and Aerospace Science and Industry (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 joint inversion model of ocean surge height and wind wave height based on the fusion of Tianmu-1 satellite-borne GNSS-R multi-modal data is shown in Table 1:

[0163]

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

[0165] Table 2 Accuracy statistics of three groups of experimental results compared with ERA5 data

[0166]

[0167] From the above table, it can be seen that the RMSE of the swell height in the case of the joint inversion model with dual output is 0.586, which is 1.72% and 9.89% lower than that in the case of the joint inversion model with single output and the DCNN model with single output, respectively. The RMSE of the wind and wave height in the case of the joint inversion model with dual output is 0.330, which is 3.28% and 0.94% lower than that in the case of the joint inversion model with single output and the DCNN model with single output, respectively. In addition, the comparison scatter density map of the swell height and wind wave height inverted by the joint inversion model with dual output and the swell height and wind wave height of ERA5, the comparison scatter density map of the swell height inverted by the joint inversion model with single output and the wind wave height inverted with single output and the swell height and wind wave height of ERA5, and the comparison scatter density map of the swell height inverted by the DCNN model with single output and the wind wave height inverted with single output and the swell height and wind wave height of ERA5 are also obtained, as shown in the figure. Figure 3-5 From the experimental results, the dual-output joint inversion model of the present invention is only slightly better than the single-output joint inversion model and the single-output DCNN model in terms of inversion accuracy, but the efficiency and economy highlight the unique advantages of this model.

[0168] Compared with many satellite-borne GNSS-R sea surface effective wave height inversion methods combined with deep learning, the advantage of the model of the present invention is that, unlike various previous models that can only invert one output target separately, this model can be directly used to output two targets at one time, and to a certain extent, the inversion accuracy of the two targets is guaranteed at the same time. This method can greatly save computing efficiency while bringing results that are not inferior to the inversion of a single output model. When applied to the economic level, it also saves costs accordingly.

[0169] In summary, the implementation of the present invention is not limited to this embodiment, but also lays a foundation for the multi-input and multi-output design of future models and opens up research ideas in related fields. This also shows that the practical significance of inverting effective wave height by satellite-borne GNSS-R technology is constantly increasing.

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

Claims

1. A joint inversion method for ocean surge height and wind wave height based on multi-modal data fusion of satellite-borne GNSS-R, characterized in that: The following steps are involved: S1, raw data acquisition, obtains the L1 data of the Tianmu-1 global navigation satellite occultation sounder-M type, the ERA5 reanalysis data set and auxiliary data as raw data; S2, preprocessing, quality control and spatiotemporal matching of raw data, and calibration of L1 data; S3, divide the data set processed in step S2 into 60% as training data, 15% as verification data, and 25% as test data; S4, construct a joint inversion model of ocean surge height and wind wave height, and train and verify the model; S5, adds an adaptive loss layer to adaptively weight the two targets of ocean surge height and wind wave height, introduces regularization on the weights to prevent overfitting and balance the output of the joint inversion model; S6, input the test data 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 wind wave height of the ERA5 data to evaluate the performance of the model.

2. The method for joint inversion of ocean surge height and wind wave height by fusion of satellite-borne GNSS-R multimodal data according to claim 1 is characterized in that: In the step S1, L1 data includes: Bistatic radar cross section, effective scattering area, raw count DDM and power DDM associated with the image; DDM related DDM bistatic radar cross section calculation coefficient, DDM peak signal-to-noise ratio, DDM power calculation coefficient, DDM waveform leading edge slope, DDM mirror reflection point normalized bistatic radar cross section, normalized mirror point signal-to-noise ratio, DDM mirror reflection point signal-to-noise ratio, UTC time corresponding to the middle moment of DDM acquisition and trajectory number of DDM waveform data; The GNSS satellite velocity X component, Y component and Z component at the signal transmission time corresponding to the DDM acquisition intermediate time associated with the transmitter and the GNSS satellite PRN code used to generate the DDM waveform; The altitude of the low-orbit satellite corresponding to the middle moment of DDM acquisition associated with the receiver, the longitude and latitude of the low-orbit satellite corresponding to the middle moment of DDM acquisition, the roll angle of the low-orbit satellite corresponding to the middle moment of DDM acquisition, and the X component, Y component, and Z component of the velocity of the low-orbit satellite corresponding to the middle moment of DDM acquisition; The height of the mirror reflection point related to the geometric relationship, the gain of the receiver antenna at the mirror reflection point, the azimuth of the mirror reflection point in the antenna coordinate system, the azimuth of the mirror reflection point in the satellite body coordinate system, the azimuth of the mirror reflection point in the orbit coordinate system, the azimuth of the mirror reflection point in the pattern coordinate system, the incident angle of the GNSS signal at the mirror reflection point, the longitude and latitude of the mirror reflection point, the altitude angle of the mirror reflection point in the antenna coordinate system, the altitude angle of the mirror reflection point in the satellite body coordinate system, the altitude angle of the mirror reflection point in the orbit coordinate system, the altitude angle of the mirror reflection point in the pattern coordinate system, and the X component, Y component and Z component of the velocity of the mirror reflection point; The ERA5 reanalysis dataset includes wind speed, wind direction, water depth, and SWH data; Auxiliary data include ERA5 precipitation data, ERA5 sea surface temperature data, SMAP sea surface salinity data, and OSCAR ocean surface current velocity data along longitude and latitude.

3. The method for joint inversion of ocean surge height and wind wave height by fusion of satellite-borne GNSS-R multimodal data according to claim 2 is characterized in that: The step S2 comprises the following sub-steps: S21, calculate characteristic observations from L1 data, perform data quality control using the L1 data product DDM quality identifier, and perform spatiotemporal matching of L1 data, ERA5 reanalysis datasets, and auxiliary data; S21, mutual calibration of characteristic observations of the four systems of L1 data: First, the raw count DDM is calibrated to the power DDM, and then the power DDM is calibrated to the bistatic radar cross section. The formula is as follows: P=C / G Where C is the raw count DDM, G is the instrument gain, and P is the power DDM. Y τ,f , are the bistatic radar cross section, power DDM, and the inverse of the effective scattering area at delay τ and Doppler f, respectively. t G t and G r is the equivalent isotropic radiated power of the GNSS transmitter and the receiver antenna gain at the specular reflection point, R t and R r are the distances from the mirror reflection point to the transmitter and receiver respectively; The characteristic observations of the four systems of L1 data are calibrated with each other using a linear model. The linear model is as follows: WITH calibrated =αZ raw +β Among them, Z raw is the uncalibrated raw observation, Z calibrated is the calibrated observation value, α is the scale factor used to adjust the scale difference between different systems, and β is the bias correction used to correct the benchmark error between different systems.

4. The method for joint inversion of ocean surge height and wind wave height by fusion of satellite-borne GNSS-R multimodal data according to claim 3 is characterized in that: In step S4, the main body of the joint inversion model combines the following five input lines: The input of input line one is the bistatic radar cross section; The input of input line 2 is GNSS-R variables and auxiliary parameters; The input of input line three is the effective scattering area; The input of input line four is the raw count DDM; The input of input line five is the power DDM.

5. The method for joint inversion of ocean surge height and wind wave height by fusion of satellite-borne GNSS-R multimodal data according to claim 4, characterized in that: The inputs of input line 1, input line 4 and input line 5 are all images, and their processing methods are as follows: 1) CNN convolution processing: Two convolutional layers are used to extract features, and then a global average pooling layer is used to reduce the dimension. The fully connected layer is then used for further feature mapping and output. This process is expressed as follows: Use two convolutional layers to extract features: X1=Conv2D(X)∈R H×W×64 X2=Conv2D(X1)∈R H×W×128 Among them, the input tensor is denoted by X∈R H×W×C Indicates that 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; Enter the global average pooling layer: X3=GlobalAveragePooling2D(X2)∈R 128 Among them, X3 is a vector of length 128, representing the extracted features, and GlobalAveragePooling2D() represents the global average pooling operation; Map X3 to 256 neurons through a fully connected layer and output X4: X4=Dense(X3)∈R 256 ; 2) Use the Transformer module to process the input data: The Transformer module uses multi-head attention and feedforward neural networks to process 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 achieve stable training and improve convergence. The feedforward network uses two fully connected layers. The first layer uses the Relu activation function to introduce nonlinearity. The second layer uses inputs.shape[-1] to project the output to the same dimension as the input. In order to prevent overfitting, random inactivation is used, and the feedforward output is added to the normalized result of the first layer to apply another residual connection. Finally, the result is normalized again as the final output of the Transformer module. The feedforward neural network here is implemented by the following formula FFN(X)=max(0,XW1+b1)W2+b2 Here W1 represents the weight matrix of the first layer; W2 represents the weight matrix of the second layer; b1, b2 represent bias terms; max(0, ) represents the Relu activation function; X represents the original input; The final output can be expressed as: Transformer(X)=LayerNorm(X+MultiHead(Q,K,V))+LayerNorm(FFN(X)) where MultiHead(Q,K,V) is a multi-head attention mechanism and LayerNorm() represents a normalization operation.

6. The method for joint inversion of ocean surge height and wind wave height by fusion of satellite-borne GNSS-R multimodal data according to claim 5, characterized in that: The input of input line three is directly processed by the Transformer module. After obtaining the final output of the Transformer module, the four lines of input line one, input line three, input line four and input line five are flattened respectively, passed through the fully connected layer respectively, and then added together to form one line.

7. The method for joint inversion of ocean surge height and wind wave height by fusion of satellite-borne GNSS-R multimodal data according to claim 6, characterized in that: The input of input line 2 uses the dynamic multi-path state space model Mamba to extract features. Three different paths are used 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 by multiplication. In addition, a gating mechanism is added, and the Softmax gating weight is used to determine the contribution of each path. The Skip connection combines the final output with the original input to promote gradient flow. In this module, assume that the input is X∈R N×D , N is the batch size, D is the feature dimension, and the final output can be expressed by the formula: Among them, Y final is the weighted sum of each path, is the result of MLP processing and one-dimensional convolution Conv1D of each path, g i represents the gating weight; After processing, it passes through a fully connected layer with 128 neurons and uses the ReLU activation function to merge with other lines. The merging process is the data fusion mechanism, which is specifically manifested as multiplication fusion. Before this step, both lines passed through a fully connected layer with 256 neurons and used the ReLU activation function to ensure the consistency of shape; After the five input lines are merged, a bidirectional convolutional long short-term memory network feature inference module is further used to effectively extract the spatial features in the input data while capturing the long-term dependencies in the time series data, thereby more comprehensively understanding the spatial-temporal characteristics of the data.

8. The method for joint inversion of ocean surge height and wind wave height by fusion of satellite-borne GNSS-R multimodal data according to claim 7, characterized in that: Step S5 includes the following sub-steps: S51, since the final output of the model is two targets, namely, swell height and wind wave height, an adaptive loss layer is added to adaptively weight the two targets. Specifically, the losses of the two targets are first calculated, and then the adaptive weights are used to update their respective weights. Finally, the total loss is selected to strengthen the focus on target one: Among them, α 1 and α 2 To adjust the weight parameters, 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 weights It is used to constrain α1 and α2 to control the complexity of the model, prevent overfitting, and balance the output of the joint inversion model. Then: The above considers both task loss and regularization loss.

9. The method for joint inversion of ocean surge height and wind wave height by fusion of satellite-borne GNSS-R multimodal data according to claim 8, characterized in that: In 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 index of high surge inversion accuracy, the relevant calculation formula is as follows: Bias=E[F(x)]-A(x) Where m is the number of data samples, y i,M is the effective surge height estimated by the model, y i,T is the swell height value from the reference dataset, y i,M and i,T The average value of E[F(x)] is the expected value of the model for a specific input x, A(x) is the true value, and A t is the true value at time t, F t is the predicted value at time t, and |·| represents an absolute value operation. The calculation method of wind and wave height assessment index is the same as above.

Images