A strong wind disaster monitoring and early warning method based on Beidou GNSS-R

Abstract

The application discloses a strong wind disaster monitoring and early warning method based on Beidou GNSS-R, comprising the following steps: 1, a UAV carrying a GNSS-R receiving device is located at the periphery of a heavy rainfall area to carry out wind speed detection operation on the water surface inside the heavy rainfall area, and data solving is performed on the obtained satellite observation; 2, an LSTM wind speed inversion model fusing double-factor decoupling is constructed, the received satellite observation and the solving result are inputted and error correction is performed, the water surface wind speed of the reflection point area is inversely predicted, and the future time step wind speed prediction is obtained based on the past time step wind speed inversion result; 3, a boundary adaptive grid and a dynamically updated Kriging interpolation model are established, the wind speed data of a missing area of a specific signal reflection point is obtained, and the wind speed data of the missing area of the reflection point is supplemented; and 4, a double-effect strong wind early warning strategy of wind speed time sequence extrapolation and inversion stability cooperation is constructed, and it is judged whether strong wind will occur and early warning is performed.

Description

Technical Field

[0001] This invention relates to the technical field of satellite-based Earth observation and remote sensing, and in particular to a method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R. Background Technology

[0002] In recent years, with the increase in outdoor activities and increasingly complex weather changes, natural disasters have become more frequent. In particular, accidents such as boat capsizing and people falling into the water due to strong winds and low visibility during activities like sightseeing, river ferries, and water sample collection are extremely difficult to resolve due to the lack of reliable meteorological information in the accident area. Given the urgent and critical situation and the safety considerations preventing personnel from venturing into dangerous areas, obtaining meteorological information about the relevant area safely, effectively, and quickly is crucial for making correct decisions and carrying out actions. Wind speed is a critical element of surface meteorology, and current observation methods mainly rely on anemometer readings and experience-based judgment based on water ripples. This is especially problematic in heavy rain scenarios, where anemometers have limited observation distances, and low clouds and rain columns further reduce visibility. Therefore, expanding the GNSS-R platform to overcome distance limitations, visual obstruction, and the effects of rainfall for surface wind speed detection has become a research hotspot in recent years.

[0003] (1) In the case of heavy rain, the wind speed is high and the visibility is low. Most of the commonly used wind measuring instruments are shore-based stations, which are only used for wind speed observation near the station location. The method of visually observing water ripples and judging wind speed by experience is also limited by the low visibility and cannot effectively obtain the water surface wind speed in the precipitation area at a distance. Sending ships or aircraft to go deep into the precipitation area to collect data is extremely dangerous, making it difficult to effectively apply the existing wind speed acquisition methods to rainstorm rescue scenarios.

[0004] (2) The location of a ship in distress is highly random and cannot be predicted in advance. Faced with complex and changing rescue needs, existing wind speed observation methods cannot achieve comprehensive and large-scale coverage, nor do they have the ability to flexibly adjust the observation area, making it difficult to achieve multi-target area and spatiotemporal continuous tracking and monitoring.

[0005] (3) For routine wind speed monitoring, setting up a large number of observation stations or increasing the number of personnel on duty is costly. In particular, when rainstorms occur, the use of contact observation may damage the instruments and equipment. Summary of the Invention

[0006] Purpose of the invention: The technical problem to be solved by the present invention is to provide a method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R, which addresses the shortcomings of the existing technology.

[0007] To address the aforementioned technical problems, this invention discloses a method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R, comprising:

[0008] Step 1: The UAV, equipped with a GNSS-R receiver, conducts wind speed detection operations on the water surface inside the heavy precipitation area from the periphery of the area. The acquired satellite observations are processed (GPS data processing is to use the signals provided by GPS satellites, receive them through a receiver, and go through a series of complex mathematical processes to obtain the position coordinates of ground points) to construct a dataset.

[0009] Step 2: Construct an LSTM wind speed inversion model that integrates two-factor decoupling. By inputting the received satellite observations and the solution results and correcting for errors, the wind speed on the water surface in the reflection point area is inverted and predicted. Based on the wind speed inversion results of the past time step, the wind speed prediction of the future time step is obtained. For example, the wind speed of the next 15 or 30 minutes can be predicted based on the wind speed of the past 30 minutes.

[0010] Step 3: Establish a boundary adaptive grid and a dynamically updated Kriging interpolation model. Use the water surface wind speed in the reflection point area for spatial adaptive interpolation to expand the inversion area and combine the wind speeds between the discrete reflection point areas, so as to obtain the wind speed data of the missing areas of specific signal reflection points and supplement the wind speed data of the missing areas of reflection points.

[0011] Step 4: Construct a dual-effect strong wind early warning strategy that combines wind speed time series extrapolation and inversion stability. The wind speed change trend is obtained by using the wind speed prediction of future time steps provided by the LSTM wind speed inversion model. Combined with the coefficient of variation of multiple continuous inverted wind speed values ​​output by the LSTM wind speed inversion model, it is determined whether strong winds will occur and an early warning is issued.

[0012] The data processing described in step 1 specifically involves:

[0013] The GNSS-R receiver analyzes the characteristic changes of L-band signals transmitted by GNSS satellites after reflection by the water surface to invert the physical parameters of the water surface. Based on the correlation between water surface roughness and wind speed, the power distribution of the reflected signal, i.e., the delayed Doppler map characteristics, is obtained by solving the problem, and the variable parameters related to water surface roughness, namely the normalized bistatic radar cross section (NBRCS) and the normalized integral delay waveform leading edge slope (LES), are obtained. These parameters, along with the reflection point location information, signal quality information, and other parameters, are used as inputs to the LSTM wind speed inversion model.

[0014] The correlation power of the delayed Doppler image obtained by solving the original satellite signal is expressed as follows:

[0015]

[0016] Where τ represents the time delay; f kThe value of the reflected signal in the selected Delayed Doppler Map (DDM) region represents the time-delayed Doppler frequency shift; f represents the Doppler frequency shift; sinc(f) is the Doppler frequency shift function; T i P represents the coherent integration time; λ represents the signal wavelength; t G represents the satellite's transmitted signal power at time t; t G is the transmit antenna gain at time t; r Indicates the receiver antenna gain; R t R is the distance between the transmitter and the mirror reflection point; r τ is the distance between the receiver and the mirror reflection point; Λ(τ) is the autocorrelation function of the pseudo-random code; σ0 is the bistatic radar cross section per unit area, i.e., the normalized bistatic radar cross section NBRCS, which reflects the sea surface scattering efficiency.

[0017] The formula for calculating the normalized bistatic radar cross section σ0 is:

[0018]

[0019] Where, p g,τ,f For the calibrated signal power, L a1 L a2 For atmospheric loss correction, l τ,f For instrument loss correction, P represents the satellite transmitted signal power. This represents the transmitter antenna gain at the reflection point. This represents the receiver antenna gain at the reflection point. For the distance loss from the transmitter to the sea surface and from the sea surface to the receiver, Λ τ；x,y S f；x,y For code correlation function and Doppler frequency shift function;

[0020] The calculation of the normalized integral delay waveform leading edge slope (LES) is as follows: Describing the steepness of the delayed Doppler waveform's leading edge, its distribution characteristics on the scattering surface represent the energy intensity distribution characteristics of the scattered signal. First, the integral delay waveform (IDW) is obtained using the delayed Doppler diagram (DDM), and the calculation expression is:

[0021]

[0022] Where N is the number of frequency units Doppler, and y is the correlation power of the sea surface scattered signal;

[0023] The normalized integral delay waveform leading edge slope (LES) is further calculated using the slope equation:

[0024]

[0025] Where, τ my represents the time delay value of the integral delay waveform IDW. m is the amplitude value of the integral delay waveform IDW, and M is the number of sampling points used to calculate the leading edge slope.

[0026] The data sources for the dataset described in step 1 include: satellite observations and meteorological element information;

[0027] The satellite observations include the longitude of the specular reflection point, the latitude of the specular reflection point, the incident angle of the specular reflection point, the distance from the transmitter to the specular reflection point, the distance from the receiver to the specular reflection point, the receiver gain at the specular reflection point, the delayed Doppler graph (DDM), the normalized bistatic radar cross section (NBRCS), and the slope of the delayed Doppler graph front (LES).

[0028] The meteorological information, used as a reference value, includes wind speed and rainfall.

[0029] Step 1 involves constructing the dataset, which includes:

[0030] Satellite observations are processed, including removing missing values, removing data from areas other than water bodies, removing data with antenna gain below a threshold of 0, and removing data with range correction gain below a threshold of 10. The formula for calculating the range correction gain (RCG) is as follows:

[0031]

[0032] in, This represents the distance from the transmitter to the point of reflection on the mirror. This represents the distance from the receiver to the point of reflection on the mirror. This represents the receiver antenna gain at the point of specular reflection.

[0033] Then, the reference wind speed and precipitation data from ERA5 (European Centre for Medium-Range Weather Forecasts, Reanalysis v5) are transformed in resolution using spatiotemporal bilinear interpolation to correspond with satellite observation data in time and space, forming a dataset for model training.

[0034] Step 2 is as follows:

[0035] Step 2-1: Use the LSTM algorithm to capture the temporal and spatial continuity of the reflection point trajectory. Combined with the fact that water surface wind speed also has continuity, train the LSTM basic wind speed inversion model. By inputting satellite observation feature parameters, obtain the preliminary predicted wind speed inversion value.

[0036] Step 2-2: Calculate the two-factor decoupling error correction for key influencing factors (distance and rainfall) in the wind speed inversion process, and obtain the distance attention correction weight and rainfall attention correction weight.

[0037] With a fixed incident angle of the satellite signal, information on signal reflection points at different horizontal distances can be obtained by changing the flight altitude of the UAV. That is, increasing the altitude allows for the detection of wind speed over water at greater distances. However, as the signal propagation path increases and is affected by rainfall, systematic biases occur in the inversion algorithm.

[0038] Steps 2-3: The fusion weight generation network takes the connection vectors of the original LSTM features, distance attention features, and rainfall attention features as input, generates normalized fusion weights through nonlinear mapping, and outputs the final wind speed prediction result.

[0039] The final wind speed inversion result is obtained by combining the basic wind speed prediction value under the LSTM algorithm framework with the two-factor error correction value. An adaptive gating mechanism is adopted to automatically learn the relative importance of distance factor and rainfall factor based on the multi-source feature state at the current moment.

[0040] The fusion weight expression is:

[0041]

[0042] in, Let t be the distance factor fusion weight. The fusion weight of the rainfall factor at time t, Let be the distance attention feature at time t. For the rainfall attention characteristics at time t, W f To fuse weights and generate the weight matrix of the network, b f The bias vector of the network is generated by fusing weights, [;;] represents the feature vector concatenation operation, H t The hidden state output obtained by the LSTM network after processing the satellite observation feature parameters at time t;

[0043] The final expression for the wind speed inversion result is:

[0044]

[0045] in, Let t be the final predicted wind speed value at time t. Linear is a linear transformation layer that maps the fused features to wind speed values.

[0046] By establishing a time-series LSTM algorithm and adding a two-factor error correction, the accuracy of wind speed inversion was improved, ensuring that the UAV platform can obtain reliable wind speed information when it is located at different observation positions. It also provides technical support for wind speed prediction in future time steps.

[0047] The two-factor decoupling error correction described in step 2-1 adopts a parallel dual-head attention architecture, with a distance attention head and a rainfall attention head designed separately;

[0048] The distance attention head takes the height sequence as its main input and maps the geometric height information into a high-dimensional feature representation, namely the distance attention feature, through a distance encoder; the distance attention calculation expression is:

[0049]

[0050] in, This is the distance attention weight matrix. These are the query matrix, key matrix, and value matrix for distance attention, respectively. These are the weight matrices for the generated query matrix, generated key matrix, and generated value matrix, respectively, representing the distance from the attention head. H t The hidden state output obtained by the LSTM network after processing the satellite observation feature parameters at time t is... Let be the distance-encoded high-dimensional feature vector at time t. For distance attention output, d k The dimension of the key vector;

[0051] The rainfall attention head takes the rainfall intensity sequence as input and adaptively adjusts the weights of each frequency component based on the rainfall intensity to obtain the rainfall attention features: The rainfall attention calculation expression is:

[0052]

[0053] in, The rainfall attention weight matrix is... These are the query matrix, key matrix, and value matrix for rainfall attention, respectively. These are the weight matrices for generating the query matrix, the key matrix, and the value matrix, respectively, for the rainfall attention head. Encode the rainfall at time t into a high-dimensional feature vector. For the output of rainfall attention features, d k is the dimension of the key vector.

[0054] Because GNSS-R technology operates in a passive satellite signal reception mode, it cannot guarantee that sufficient signal reflection points are always located in a specific observation area. Therefore, a spatiotemporal interpolation model is needed to supplement wind speed data in areas lacking signal reflection points. Utilizing the maneuverability of UAV platforms, rapid detection of different bodies of water can be achieved. By changing observation positions, as many regional surface wind speed inversion results as possible can be obtained, supplementing wind speed data in areas lacking reflection points. This "missing" refers to areas where no satellite signal reflection point exists at any given moment, making it impossible to directly obtain wind speed data through "signal inversion." For example, in a mission targeting a lake, if at a certain moment only the southwest corner has a reflection point, then the northeast corner of the lake is the "missing reflection point area." This supplementation is achieved through steps 3-1 and 3-2, which can be simply understood as "data interpolation under certain rules" (inserting new data into areas lacking data based on the characteristics of existing data; the new data does not depend on reflection point inversion).

[0055] By establishing a boundary adaptive grid kriging interpolation model based on classical kriging interpolation theory, considering the irregular characteristics of water body boundaries, the physical characteristics of wind fields, and the confidence and high-frequency update characteristics oriented towards rescue needs, a boundary adaptive grid partitioning mechanism and physical constraints are introduced, and uncertainty quantification and adaptive update blocks are added to achieve the best spatial coverage effect with the least observation cost.

[0056] Step 3 specifically involves:

[0057] Step 3-1: Based on the classic Kriging interpolation model, introduce physical constraints on wind speed variation to construct a semi-variogram model that conforms to the spatial variation law of wind speed.

[0058] Step 3-2: Construct a boundary adaptive grid system. The grid serves as the basic spatial unit for wind speed interpolation. Each grid contains a specific wind speed value, representing the wind speed within that grid region, expressed as:

[0059] G a ={T k |T k =Delaunay(P b ∪P o ,ρ k )}

[0060] Among them, G a For an adaptive grid system, T k For the k-th triangular mesh element, P b Let P be the set of points constrained by the water boundary. o Let ρ be the set of observation point locations. k Let be the grid density parameter for the k-th region;

[0061] Further, an adaptive grid density adjustment function is established to adjust the grid density and size based on the spatial distribution relationship between the region with known wind speed and the region where the wind speed needs to be interpolated, and in conjunction with the shape of the water body boundary:

[0062] ρ(s)=ρ0·exp(-α·d o (s))·(1+β c ·C b (s))

[0063] Where ρ(s) is the grid density at position s, ρ0 is the base grid density, α is the observation point distance attenuation coefficient, and d o (s) is the distance from position s to the nearest observation point, β c For boundary complexity weights, C b (s) is the boundary complexity index at position s;

[0064] Step 3-3: Based on the semi-variogram model and boundary adaptive grid system that conform to the spatial variation law of wind speed, construct an adaptive grid Kriging wind speed interpolation model, the mathematical expression of which is:

[0065]

[0066] in, V(s) is the interpolated predicted wind speed at location s0. i ) represents the i-th observation point s i The inverted wind speed value at λ i ρ(s) represents the Kriging weight coefficient for the i-th observation point, n is the total number of observation points participating in the interpolation, and ρ(s) represents the weight coefficient for the i-th observation point. i ) represents the position s i The grid density at that location, where μ is the Lagrange multiplier;

[0067] The wind speed interpolation prediction result is obtained based on the adaptive grid Kriging wind speed interpolation model, and the uncertainty estimate is provided through the prediction variance. The expression for calculating the prediction variance is as follows:

[0068]

[0069] in, Let be the Kriging prediction variance at position s0, and γ(0) be the semivariogram value at zero distance;

[0070] Steps 3-4: After realizing the spatial interpolation of wind speed, the Kriging model can provide the prediction variance of each interpolation point. Based on the prediction variance, a confidence interval is constructed and the confidence level is calculated. When the prediction uncertainty at the grid point is large, it indicates that the area needs to be supplemented with observations, providing a reference for the UAV to change to a new observation point.

[0071] The expression for calculating grid confidence intervals is:

[0072]

[0073] Among them, CI α (s0) represents the confidence interval for position s with a confidence level of (1-α). Let z be the Kriging prediction value at position s0. α / 2 σ is the two-sided critical value of the standard normal distribution. k (s0) represents the standard deviation of the Kriging prediction at location s0;

[0074] Steps 3-5: After obtaining the adaptive grid kriging wind speed interpolation results and confidence levels, perform real-time weight updates and adjustments, expressed as follows:

[0075]

[0076] in, The weight coefficients are updated at time t+1. Let L be the weight coefficient at time t, η be the learning rate parameter, and L be the weight coefficient at time t. t This is the total loss function that includes prediction error and physical constraints;

[0077] Furthermore, spatial correlation adaptive updates are performed to obtain a semi-variogram function that changes synchronously with the newly retrieved wind speed data, enabling the grid interpolated wind speed to maintain a high-frequency update frequency. The expression is:

[0078] γ u (h)=(1-ω)·γ o (h)+ω·γ e (h)

[0079] Where, γ u (h) is the updated semi-mutation function, γ o (h) is the semi-variogram function before the update, γ e (h) is the empirical semivariogram based on new observation data, and ω is the update weight used to control the fusion ratio of new and old information, with a value range of 0 to 1.

[0080] Traditional Kriging interpolation algorithms can effectively capture the spatial correlation of meteorological elements. To address the actual needs of UAV platforms in retrieving wind speeds over inland waters, we have improved the interpolation grid and added physical constraints, or filled in the data with actual wind speed data to achieve wind speed information acquisition for multiple target areas. At the same time, we have added a module for confidence assessment and adaptive update of interpolation results to provide rescue personnel with more realistic decision-making assistance.

[0081] Step 3-1 is as follows:

[0082] Step 3-1-1: The Kriging interpolation model for the boundary adaptive mesh is a targeted improvement on the classic Kriging interpolation model. Its basic expression is:

[0083]

[0084] in, V(s) is the interpolated predicted wind speed at location s0. i ) represents the i-th observation point s i The measured wind speed value at the location, λ i is the Kriging weight coefficient for the i-th observation point, n is the total number of observation points participating in the interpolation, is the Lagrange multiplier, and μ is used to satisfy the unbiased estimation condition;

[0085] The systematic expression for weight calculation in the classic Kriging interpolation model is:

[0086] K·λ=k0

[0087] Expanded to:

[0088]

[0089] Where K is an (n+1)×(n+1) dimensional semivariogram matrix, λ is the weight coefficient vector, k0 is the target vector, and γ(s) i ,s0) is the position s i The semivariogram values ​​between s0 and s0;

[0090] Step 3-1-2: Based on the classic Kriging interpolation model, physical constraints on wind speed variation are introduced to obtain a semi-variogram model that conforms to the spatial variation law of wind speed, expressed as:

[0091]

[0092] Where γ(h) is the semivariogram value at a distance h, C0 is the nugget effect, C1 is the partial sill value, a is the range parameter, h is the spatial distance between observation points, β is the physical constraint weighting coefficient, and Φ p (h) is a constraint function based on physical laws.

[0093] Step 4 specifically involves:

[0094] Faced with rapidly changing wind speeds and the possibility of strong winds, and given the availability of spatiotemporally continuous wind speed information, forward prediction is made by analyzing past time step wind speed changes to perceive future wind speeds and their trends. Combining the inherent characteristic of decreased accuracy as wind speed increases during GNSS-R high-wind-speed inversion, a sliding window is used to calculate the coefficient of variation of the inversion results. This constructs a dual-effect strong wind early warning strategy that coordinates wind speed time series extrapolation and inversion stability, providing effective early warning of strong wind disasters.

[0095] Step 4-1: Obtain the time-series extrapolated wind speed (i.e., the wind speed at future times) through LSTM multi-step forward prediction. The expression is:

[0096]

[0097] in, Let h be the predicted wind speed at time t+k. t Let N be the LSTM hidden state at time t, and N be the prediction time step.

[0098] Step 4-2: Based on the predicted wind speed, calculate the rate of change of wind speed, expressed as:

[0099]

[0100] Δ 2 V t+k =ΔV t+k -ΔV t+k-1

[0101] Where, ΔV t+k Let Δ be the rate of change of wind speed at time t+k. 2 V t+k Let be the wind speed acceleration at time t+k;

[0102] Step 4-3: Record the wind speed values ​​output by the wind speed inversion model from Step 2-3 at multiple consecutive time points. Use the sliding window coefficient of variation method to determine the stability of the wind speed inversion. Calculate the wind speed coefficient of variation within the time window based on the wind speed values ​​at multiple consecutive time points. A large coefficient of variation indicates rapid wind speed changes, while a small coefficient indicates relatively stable wind speed with smaller fluctuations. The expression for calculating the sliding window coefficient of variation is:

[0103]

[0104] Among them, CV t Let σ be the coefficient of variation at time t. t v represents the standard deviation of wind speed within the sliding window. t V is the average wind speed within the sliding window, W is the length of the sliding window, and V is the average wind speed within the sliding window. i The actual wind speed inversion value at time i;

[0105] Step 4-4: Combine the wind speed and wind speed change rate predicted by wind speed time series extrapolation with the coefficient of variation of the sliding time window to conduct a coordinated strong wind warning. In actual use, the standard should be formulated and implemented in accordance with relevant documents or regulations.

[0106] The definitions of the proper nouns used in this invention include:

[0107] GNSS (Global Navigation Satellite System)

[0108] GNSS-R Global Navigation Satellite System Reflectometry

[0109] Delay-Doppler Map (DDM)

[0110] NBRCS (Normalized Bistatic Radar Cross Section)

[0111] LES Normalized Integral Delay Waveform Leading Edge Slope

[0112] LSTM (Long Short-Term Memory) network

[0113] IDW Integral Delay Waveform

[0114] RCG Range Corrected Gain

[0115] ERA5 (European Centre for Medium-Range Weather Forecasts) fifth-generation reanalysis data (ECMWF Reanalysis v5)

[0116] Beneficial effects:

[0117] Technical aspects:

[0118] 1. Breaking through traditional monitoring limitations: By using an innovative method of GNSS-R receiving equipment on a drone, the technical difficulties of traditional wind measuring instruments, such as limited observation distance, visual obstruction, and the impact of rainfall, have been overcome, enabling multi-directional, long-distance, and all-weather water surface wind speed detection.

[0119] 2. High-precision wind speed inversion technology: An LSTM wind speed inversion model with dual-factor decoupling was constructed. Through dual error correction of distance attention and rainfall attention, the accuracy and stability of wind speed inversion under heavy precipitation observation conditions were significantly improved.

[0120] 3. Enhanced spatial coverage: An adaptive boundary grid and dynamically updated Kriging interpolation model were established to supplement wind speed data in areas where signal reflection points are missing, thereby improving the utilization rate of satellite observation data, expanding the monitoring coverage, and providing spatially continuous wind speed information.

[0121] 4. Intelligent early warning strategy: A dual-effect strong wind early warning strategy combining wind speed time series extrapolation and inversion stability was constructed. By combining future wind speed prediction and coefficient of variation analysis, a more accurate and timely strong wind disaster early warning was achieved.

[0122] Application level:

[0123] 1. Addressing emergency rescue needs: In response to emergencies such as ships in distress and people falling into the water during heavy rain, it provides a safe and rapid means of obtaining meteorological information of the accident area, providing key support for rescue decision-making.

[0124] 2. Enhance monitoring flexibility: It has the ability to flexibly adjust the observation area and can conduct multi-target area and spatiotemporal continuous tracking and monitoring according to actual needs, adapting to complex and ever-changing rescue needs.

[0125] 3. Reduce monitoring costs and risks: It avoids the high cost of setting up a large number of observation stations or increasing the number of on-duty personnel, reduces the risk of damage to instruments and equipment during rainstorm disasters, and realizes non-contact safety monitoring.

[0126] 4. Wide range of applications: It can be applied to accident rescue such as falling into water and capsizing of boats, as well as to the safety of various water-related activities such as river ferries and scenic sightseeing, and has good value for promotion and application. Attached Figure Description

[0127] Figure 1 This is a flowchart of the algorithm of the present invention. Detailed Implementation

[0128] A method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R, comprising:

[0129] Step 1: The UAV, equipped with a GNSS-R receiver, is located outside the area of ​​heavy precipitation to conduct wind speed detection on the water surface inside the area. The acquired satellite observations are processed to construct a dataset.

[0130] Step 2: Construct an LSTM wind speed inversion model that integrates two-factor decoupling. By inputting the received satellite observations and the solution results and correcting for errors, the wind speed at the water surface in the reflection point area is inverted and predicted. Based on the wind speed inversion results of the past time step, the wind speed prediction for the future time step is obtained.

[0131] Step 3: Establish a boundary adaptive grid and a dynamically updated Kriging interpolation model. Use the water surface wind speed in the reflection point area for spatial adaptive interpolation to expand the inversion area and combine the wind speeds between the discrete reflection point areas, so as to obtain the wind speed data of the missing areas of specific signal reflection points and supplement the wind speed data of the missing areas of reflection points.

[0132] Step 4: Construct a dual-effect strong wind early warning strategy that combines wind speed time series extrapolation and inversion stability. The wind speed change trend is obtained by using the wind speed prediction of future time steps provided by the LSTM wind speed inversion model. Combined with the coefficient of variation of multiple continuous inverted wind speed values ​​output by the LSTM wind speed inversion model, it is determined whether strong winds will occur and an early warning is issued.

[0133] The data processing described in step 1 specifically involves:

[0134] The GNSS-R receiver analyzes the characteristic changes of L-band signals transmitted by GNSS satellites after reflection by the water surface to invert the physical parameters of the water surface. Based on the correlation between water surface roughness and wind speed, it calculates the power distribution of the reflected signal, i.e., the delayed Doppler map characteristics, and obtains the variable parameters related to water surface roughness: the normalized bistatic radar cross section (NBRCS) and the normalized integral delay waveform leading edge slope (LES), and then combines them with the reflection point location information and signal quality information.

[0135] The correlation power of the delayed Doppler image obtained by solving the original satellite signal is expressed as follows:

[0136]

[0137] Where τ represents the time delay; f k The value of the reflected signal in the selected time-delay Doppler shift region corresponds to the time-delay Doppler map region; f represents the Doppler frequency shift; sinc(f) is the Doppler frequency shift function; T i P represents the coherent integration time; λ represents the signal wavelength; t G represents the satellite's transmitted signal power at time t; t G is the transmit antenna gain at time t; r Indicates the receiver antenna gain; R t R is the distance between the transmitter and the mirror reflection point; r τ is the distance between the receiver and the mirror reflection point; Λ(τ) is the autocorrelation function of the pseudo-random code; σ0 is the bistatic radar cross section per unit area, i.e., the normalized bistatic radar cross section NBRCS.

[0138] The formula for calculating the normalized bistatic radar cross section σ0 is:

[0139]

[0140] Where, p g,τ,f For the calibrated signal power, L a1 L a2 For atmospheric loss correction, l τ,f For instrument loss correction, P represents the satellite transmitted signal power. This represents the transmitter antenna gain at the reflection point. This represents the receiver antenna gain at the reflection point. For the distance loss from the transmitter to the sea surface and from the sea surface to the receiver, Λ τ；x,y S f；x,y For code correlation function and Doppler frequency shift function;

[0141] The calculation of the normalized integral delay waveform leading edge slope (LES) is as follows: First, the integral delay waveform (IDW) is obtained using the delayed Doppler plot, and the calculation expression is:

[0142]

[0143] Where N is the number of frequency units Doppler, and y is the correlation power of the sea surface scattered signal;

[0144] The normalized integral delay waveform leading edge slope (LES) is further calculated using the slope equation:

[0145]

[0146] Where, τ m y represents the time delay value of the integral delay waveform IDW. m is the amplitude value of the integral delay waveform IDW, and M is the number of sampling points used to calculate the leading edge slope.

[0147] The data sources for the dataset described in step 1 include: satellite observations and meteorological element information;

[0148] The satellite observations include the longitude of the specular reflection point, the latitude of the specular reflection point, the incident angle of the specular reflection point, the distance from the transmitter to the specular reflection point, the distance from the receiver to the specular reflection point, the receiver gain at the specular reflection point, the delayed Doppler plot, the normalized bistatic radar cross section, and the slope of the leading edge of the delayed Doppler plot.

[0149] The meteorological information includes wind speed and rainfall.

[0150] Step 1 involves constructing the dataset, which includes:

[0151] Satellite observations are processed, including removing missing values, removing data from areas other than water bodies, removing data with antenna gain below a threshold, and removing data with range correction gain below a threshold. The formula for calculating the range correction gain (RCG) is as follows:

[0152]

[0153] in, This represents the distance from the transmitter to the point of reflection on the mirror. This represents the distance from the receiver to the point of reflection on the mirror. This represents the receiver antenna gain at the point of specular reflection.

[0154] Then, the reference wind speed and rainfall data from ERA5 are transformed in resolution using spatiotemporal bilinear interpolation to correspond with satellite observation data in time and space, forming a dataset for model training.

[0155] Step 2 is as follows:

[0156] Step 2-1: Use the LSTM algorithm to capture the temporal and spatial continuity of the reflection point trajectory. Combined with the fact that water surface wind speed also has continuity, train the LSTM basic wind speed inversion model. By inputting satellite observation feature parameters, obtain the preliminary predicted wind speed inversion value.

[0157] Step 2-2: Calculate the two-factor decoupling error correction for key influencing factors in the wind speed inversion process, and obtain the distance attention correction weight and rainfall attention correction weight.

[0158] Steps 2-3: The fusion weight generation network takes the connection vectors of the original LSTM features, distance attention features, and rainfall attention features as input, generates normalized fusion weights through nonlinear mapping, and outputs the final wind speed prediction result.

[0159] The fusion weight expression is:

[0160]

[0161] in, Let t be the distance factor fusion weight. The fusion weight of the rainfall factor at time t, Let be the distance attention feature at time t. For the rainfall attention characteristics at time t, W f To fuse weights and generate the weight matrix of the network, b f The bias vector of the network is generated by fusing weights, [;;] represents the feature vector concatenation operation, H t The hidden state output obtained by the LSTM network after processing the satellite observation feature parameters at time t;

[0162] The final expression for the wind speed inversion result is:

[0163]

[0164] in, Let t be the final predicted wind speed value at time t. Linear is a linear transformation layer that maps the fused features to wind speed values.

[0165] The two-factor decoupling error correction described in step 2-1 adopts a parallel dual-head attention architecture, with a distance attention head and a rainfall attention head designed separately;

[0166] The distance attention head takes the height sequence as its main input and maps the geometric height information into a high-dimensional feature representation, namely the distance attention feature, through a distance encoder; the distance attention calculation expression is:

[0167]

[0168] in, This is the distance attention weight matrix. These are the query, key, and value matrices for distance attention. These are the weight matrices for distance from the attention head, H. t The hidden state output obtained by the LSTM network after processing the satellite observation feature parameters at time t is... Let be the distance-encoded high-dimensional feature vector at time t. For distance attention output, d k The dimension of the key vector;

[0169] The rainfall attention head takes the rainfall intensity sequence as input and adaptively adjusts the weights of each frequency component based on the rainfall intensity to obtain the rainfall attention features: The rainfall attention calculation expression is:

[0170]

[0171] in, The rainfall attention weight matrix is... These are the query, key, and value matrices for rainfall attention. These are the weight matrices for the rainfall attention head. Encode the rainfall at time t into a high-dimensional feature vector. For the output of rainfall attention features, d k is the dimension of the key vector.

[0172] Step 3 specifically involves:

[0173] Step 3-1: Based on the classic Kriging interpolation model, introduce physical constraints on wind speed variation to construct a semi-variogram model that conforms to the spatial variation law of wind speed.

[0174] Step 3-2: Construct a boundary adaptive grid system. The grid serves as the basic spatial unit for wind speed interpolation. Each grid contains a specific wind speed value, representing the wind speed within that grid region, expressed as:

[0175] G a ={T k |T k =Delaunay(P b ∪P o ,ρ k)}

[0176] Among them, G a For an adaptive grid system, T k For the k-th triangular mesh element, P b Let P be the set of points constrained by the water boundary. o Let ρ be the set of observation point locations. k Let be the grid density parameter for the k-th region;

[0177] Further, an adaptive grid density adjustment function is established to adjust the grid density and size based on the spatial distribution relationship between the region with known wind speed and the region where the wind speed needs to be interpolated, and in conjunction with the shape of the water body boundary:

[0178] ρ(s)=ρ0·exp(-α·d o (s))·(1+β c ·C b (s))

[0179] Where ρ(s) is the grid density at position s, ρ0 is the base grid density, α is the observation point distance attenuation coefficient, and d o (s) is the distance from position s to the nearest observation point, β c For boundary complexity weights, C b (s) is the boundary complexity index at position s;

[0180] Step 3-3: Based on the semi-variogram model and boundary adaptive grid system that conform to the spatial variation law of wind speed, construct an adaptive grid Kriging wind speed interpolation model, the mathematical expression of which is:

[0181]

[0182] in, V(s) is the interpolated predicted wind speed at location s0. i ) represents the i-th observation point s i The inverted wind speed value at λ i ρ(s) represents the Kriging weight coefficient for the i-th observation point, n is the total number of observation points participating in the interpolation, and ρ(s) represents the weight coefficient for the i-th observation point. i ) represents the position s i The grid density at that location, where ρ is the Lagrange multiplier;

[0183] The wind speed interpolation prediction result is obtained based on the adaptive grid Kriging wind speed interpolation model, and the uncertainty estimate is provided through the prediction variance. The expression for calculating the prediction variance is as follows:

[0184]

[0185] in, Let be the Kriging prediction variance at position s0, and γ(0) be the semivariogram value at zero distance;

[0186] Steps 3-4: After realizing the spatial interpolation of wind speed, the Kriging model can provide the prediction variance of each interpolation point. Based on the prediction variance, a confidence interval is constructed and the confidence level is calculated. When the prediction uncertainty at the grid point is large, it indicates that the area needs to be supplemented with observations, providing a reference for the UAV to change to a new observation point.

[0187] The expression for calculating grid confidence intervals is:

[0188]

[0189] Among them, CI α (s0) represents the confidence interval for position s with a confidence level of (1-α). Let z be the Kriging prediction value at position s0. α / 2 σ is the two-sided critical value of the standard normal distribution. k (s0) represents the standard deviation of the Kriging prediction at location s0;

[0190] Steps 3-5: After obtaining the adaptive grid kriging wind speed interpolation results and confidence levels, perform real-time weight updates and adjustments, expressed as follows:

[0191]

[0192] in, The weight coefficients are updated at time t+1. Let L be the weight coefficient at time t, η be the learning rate parameter, and L be the weight coefficient at time t. t This is the total loss function that includes prediction error and physical constraints;

[0193] Furthermore, spatial correlation adaptive updates are performed to obtain a semi-variogram function that changes synchronously with the newly retrieved wind speed data, enabling the grid interpolated wind speed to maintain a high-frequency update frequency. The expression is:

[0194] γ u (h)=(1-ω)·γ o (h)+ω·γ e (h)

[0195] Where, γ u (h) is the updated semi-mutation function, γ o (h) is the semi-variogram function before the update, γ e (h) is the empirical semivariogram based on new observation data, and ω is the update weight used to control the fusion ratio of new and old information, with a value range of 0 to 1.

[0196] Step 3-1 is as follows:

[0197] Step 3-1-1: The Kriging interpolation model for the boundary adaptive mesh is a targeted improvement on the classic Kriging interpolation model. Its basic expression is:

[0198]

[0199] in, V(s) is the interpolated predicted wind speed at location s0. i ) represents the i-th observation point s i The measured wind speed value at the location, λ i is the Kriging weight coefficient for the i-th observation point, n is the total number of observation points participating in the interpolation, is the Lagrange multiplier, and μ is used to satisfy the unbiased estimation condition;

[0200] The systematic expression for weight calculation in the classic Kriging interpolation model is:

[0201] K·λ=k0

[0202] Expanded to:

[0203]

[0204] Where K is an (n+1)×(n+1) dimensional semivariogram matrix, λ is the weight coefficient vector, k0 is the target vector, and γ(s) i ,s0) is the position s i The semivariogram values ​​between s0 and s0;

[0205] Step 3-1-2: Based on the classic Kriging interpolation model, physical constraints on wind speed variation are introduced to obtain a semi-variogram model that conforms to the spatial variation law of wind speed, expressed as:

[0206]

[0207] Where γ(h) is the semivariogram value at a distance h, C0 is the nugget effect, C1 is the partial sill value, a is the range parameter, h is the spatial distance between observation points, β is the physical constraint weighting coefficient, and Φ p (h) is a constraint function based on physical laws.

[0208] Step 4 specifically involves:

[0209] Step 4-1: Obtain the time-series extrapolated wind speed through LSTM multi-step forward prediction, expressed as:

[0210]

[0211] in, Let h be the predicted wind speed at time t+k. tLet N be the LSTM hidden state at time t, and N be the prediction time step.

[0212] Step 4-2: Based on the predicted wind speed, calculate the rate of change of wind speed, expressed as:

[0213]

[0214] Δ 2 V t+k =ΔV t+k -ΔV t+k-1

[0215] Where, ΔV t+k Let Δ be the rate of change of wind speed at time t+k. 2 V t+k Let be the wind speed acceleration at time t+k;

[0216] Step 4-3: Record the wind speed values ​​output by the wind speed inversion model from Step 2-3 at consecutive time points. Use the sliding window coefficient of variation method to determine the stability of the wind speed inversion. Calculate the wind speed coefficient of variation within the time window based on the wind speed values ​​at multiple consecutive time points. A large coefficient of variation indicates rapid wind speed changes, while a small coefficient indicates relatively stable wind speed with smaller fluctuations. The expression for calculating the sliding window coefficient of variation is:

[0217]

[0218] Among them, CV t Let σ be the coefficient of variation at time t. t Let μ be the standard deviation of wind speed within the sliding window. t V is the average wind speed within the sliding window, W is the length of the sliding window, and V is the average wind speed within the sliding window. i The actual wind speed inversion value at time i;

[0219] Step 4-4: Combine the wind speed and wind speed change rate predicted by wind speed time series extrapolation with the coefficient of variation of the sliding time window to conduct a coordinated strong wind warning.

[0220] Example 1:

[0221] This embodiment uses the strong wind early warning and monitoring of sightseeing boats on a lake in a natural scenic area as an example to explain in detail the implementation process of the present invention:

[0222] Step 1: Acquire and process satellite observation data to obtain the feature inputs for the wind speed inversion model:

[0223] A multi-rotor UAV equipped with a GNSS-R receiver was selected, with its flight altitude set at 200 meters, and positioned in a safe area outside the rain belt of the target lake. The target monitoring area was an elliptical lake with an area of ​​approximately 15 square kilometers, in which three tourist sightseeing boats were operating.

[0224] The GNSS-R receiver receives direct and reflected satellite signals in the area, obtains the delayed Doppler map through raw data processing, and extracts key parameters, including the coordinates of the specular reflection point, the incident angle, the normalized bistatic radar cross section (NBRCS) value, the normalized integral delay waveform leading edge slope (LES) value, and the distances from the reflection point to the signal transmitter and receiver, respectively.

[0225] Step 2: Based on satellite observation data, perform LSTM wind speed inversion using a fusion of two factors and decoupled methods to obtain the inverted wind speed at the signal reflection point:

[0226] An LSTM network with 64 hidden units was constructed, and the time series length was set to 10 time steps. The encoding dimension of the distance attention head was set to 32, and the encoding dimension of the rainfall attention head was set to 16. The fusion weight generation network adopted a 2-layer fully connected structure with a learning rate of 0.001.

[0227] Under the current observation conditions (drone altitude 200m, horizontal distance approximately 280m, rainfall intensity 20mm / h), the model outputs a preliminary wind speed inversion value of 9.2m / s. After two-factor error correction, with a distance factor weight of 0.45 and a rainfall factor weight of 0.55, the final wind speed inversion result is 8.4m / s.

[0228] Step 3: Divide the lake surface into boundary adaptive grids and perform dynamically updated improved Kriging model wind speed interpolation to expand the lake surface wind speed monitoring range and improve the spatiotemporal resolution of monitoring.

[0229] Based on the irregular boundary characteristics of the lake area, an adaptive mesh system containing 458 triangular mesh elements was constructed. The basic mesh density ρ0 was set to 0.025, the observation point distance attenuation coefficient a was set to 0.18, and the boundary complexity weight β was set to... c Set it to 0.35.

[0230] Since the northeastern corner of the lake area lacks a signal reflection point at the current moment, and the reflection is mainly concentrated in the central part of the lake area, an improved Kriging interpolation model is used to estimate and assign grid wind speeds in this area. The interpolation results show that the wind speed in this area is 8.1 m / s, with a confidence interval of [7.6, 8.6] m / s and a confidence level of 92%.

[0231] Step 4: Issue a strong wind warning based on changes in wind speed:

[0232] A sliding window of 10 minutes was set, and an LSTM time series model was used to predict the wind speed trend over the next 30 minutes. The prediction results showed that the wind speed would reach 12.1 m / s after 15 minutes and 15.7 m / s after 30 minutes, with a wind speed change rate of 0.24 m / s. 2.

[0233] Meanwhile, the coefficient of variation (CV) of wind speed within the current 10-minute window is 0.12, exceeding the set threshold of 0.1, indicating that the wind speed is changing rapidly. Based on the combined time-series extrapolation results and CV analysis, the system predicts that strong winds (≥12 m / s) will occur within 30 minutes. Furthermore, considering the geographical environment of the scenic area and the wave resistance of the sightseeing boats, warning levels are established, including Level 1 (wind speed ≥18 m / s, extremely high probability of boat capsizing, immediate rendezvous with the nearest shore or return to the dock), Level 2 (wind speed ≥15 m / s, relatively high probability of boat capsizing, immediate orderly return to the dock), and Level 3 (wind speed ≥12 m / s, risk of boat capsizing, limiting the distance of boats from shore). The drone automatically sends a Level 2 warning signal to the scenic area management duty room based on the wind speed inversion monitoring results.

[0234] Warning execution result:

[0235] The system simultaneously sent an early warning to the lake management center and three tourist boats: "The wind speed in the lake area is expected to reach 15.7 m / s within 30 minutes. The boats are required to immediately organize their formation and return to port to avoid the wind." After receiving the warning, the boats returned to port in time, successfully avoiding potential safety accidents caused by the strong winds.

[0236] Actual verification shows that the wind speed prediction error of this warning is less than 1.2 m / s, and the warning time lead time reaches 28 minutes, providing sufficient time for the safe evacuation of ships.

[0237] This invention provides a method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R. Many methods and approaches exist for implementing this technical solution; the above description is merely a preferred embodiment. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of this invention, and these improvements and modifications should also be considered within the scope of protection of this invention. All components not explicitly stated in this embodiment can be implemented using existing technologies.

Claims

1. A method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R, characterized in that, include: Step 1: The UAV, equipped with a GNSS-R receiver, is located outside the area of ​​heavy precipitation to conduct wind speed detection on the water surface inside the area. The acquired satellite observations are processed to construct a dataset. Step 2: Construct a LSTM wind speed inversion model that integrates two-factor decoupling. By inputting the received satellite observations and the solution results and applying error correction, the model inverts and predicts the surface wind speed in the reflection point area. Based on the wind speed inversion results of past time steps, the model predicts the wind speed for future time steps, including: Step 2-1: Use the LSTM algorithm to capture the temporal and spatial continuity of the reflection point trajectory. Combined with the fact that water surface wind speed also has continuity, train the LSTM basic wind speed inversion model. By inputting satellite observation feature parameters, obtain the preliminary predicted wind speed inversion value. Step 2-2 involves calculating the two-factor decoupling error correction for key influencing factors in the wind speed inversion process, and obtaining the distance attention correction weight and rainfall attention correction weight. These key influencing factors include distance and rainfall. Steps 2-3: The fusion weight generation network takes the connection vectors of the original LSTM features, distance attention features, and rainfall attention features as input, generates normalized fusion weights through nonlinear mapping, and outputs the final wind speed prediction result. The two-factor decoupling error correction described in step 2-2 adopts a parallel dual-head attention architecture, with a distance attention head and a rainfall attention head designed separately; Step 3 involves establishing a boundary adaptive grid and a dynamically updated Kriging interpolation model. Spatial adaptive interpolation is performed using the surface wind speed in the reflection point region to expand the inversion area and combine wind speeds between discrete reflection point regions. This obtains wind speed data for areas lacking specific signal reflection points, supplementing the missing wind speed data in these areas. Specifically, Step 3 includes: Step 3-1: Based on the classic Kriging interpolation model, introduce physical constraints on wind speed variation to construct a semi-variogram model that conforms to the spatial variation law of wind speed. Step 3-2: Construct a boundary adaptive grid system. The grid serves as the basic spatial unit for wind speed interpolation. Each grid contains a specific wind speed value, representing the wind speed within that grid region, expressed as: in, For an adaptive grid system, For the first A triangular grid cell, For the set of constraint points of the water boundary, For the set of observation point locations, For the first The grid density parameter for each region, Delaunay represents the construction of triangular grid cells; An adaptive grid density adjustment function is established to adjust the grid density and size based on the spatial distribution relationship between the region with known wind speed and the region where the wind speed needs to be interpolated, and in conjunction with the shape of the water body boundary. in, For position Mesh density at that location, Based on grid density, The distance attenuation coefficient is the distance between the observation points. For position Distance to the nearest observation point For boundary complexity weights, For position The boundary complexity index at the location; Step 3-3: Based on the semi-variogram model and boundary adaptive grid system that conform to the spatial variation law of wind speed, construct an adaptive grid Kriging wind speed interpolation model, the mathematical expression of which is: in, For the location The interpolated predicted wind speed value, For the first Observation points The inversion wind speed value at that location, For the first Kriging weight coefficients for each observation point The total number of observation points participating in the interpolation. For position Mesh density at that location, For position Mesh density at that location, It is a Lagrange multiplier; The wind speed interpolation prediction result is obtained based on the adaptive grid Kriging wind speed interpolation model, and the uncertainty estimate is provided through the prediction variance. The expression for calculating the prediction variance is as follows: in, For position Kriging's prediction variance The value of the semivariogram at the zero distance; Steps 3-4: After performing spatial interpolation of wind speed, the Kriging model is used to provide the prediction variance for each interpolation point. Based on the prediction variance, grid confidence intervals are constructed to calculate the confidence level. The expression for calculating grid confidence intervals is: in, For position The level of handling information is The confidence interval, For position Kriging's predicted value, This is the two-sided critical value of the standard normal distribution. For position Kriging's standard deviation of prediction; Steps 3-5: After obtaining the adaptive grid kriging wind speed interpolation results and confidence levels, perform real-time weight updates and adjustments, expressed as follows: in, For the first The weighting coefficients are updated in real time. For the first Weighting coefficients at time points, This is the learning rate parameter, with a value ranging from 0.01 to 0.

05. This is the total loss function that includes prediction error and physical constraints; Furthermore, spatial correlation adaptive updates are performed to obtain a semi-variogram function that changes synchronously with the newly retrieved wind speed data, enabling the grid interpolated wind speed to maintain a high frequency of updates. The expression is: in, For the updated semi-mutation function, The semi-mutation function before the update. For the empirical semivariogram based on new observation data, The updated weights are used to control the fusion ratio of new and old information, and their values ​​range from 0 to 1. Step 4: Construct a dual-effect strong wind early warning strategy that combines wind speed time series extrapolation and inversion stability. The wind speed change trend is obtained by using the wind speed prediction of future time steps provided by the LSTM wind speed inversion model. Combined with the coefficient of variation of multiple continuous inverted wind speed values ​​output by the LSTM wind speed inversion model, it is determined whether strong winds will occur and an early warning is issued.

2. The method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R according to claim 1, characterized in that, The data processing described in step 1 specifically involves: The GNSS-R receiver analyzes the characteristic changes of L-band signals transmitted by GNSS satellites after reflection by the water surface to invert the physical parameters of the water surface. Based on the correlation between water surface roughness and wind speed, the power distribution of the reflected signal, i.e., the delayed Doppler map characteristics, is obtained by solving the problem, and the variable parameters related to water surface roughness are obtained: normalized bistatic radar cross section (NBRCS) and normalized integral delay waveform leading edge slope (LES). These parameters, along with the reflection point location information and signal quality information, are used as inputs to the LSTM wind speed inversion model.

3. The method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R according to claim 2, characterized in that, The correlation power of the delayed Doppler image is expressed as follows: in, Indicates time delay; This represents the reflected signal power value of the selected time-delay Doppler frequency shift corresponding to the delay Doppler map region; Indicates Doppler frequency shift; This is the Doppler frequency shift function; Indicates the coherent integration time; Indicates the signal wavelength; This represents the satellite's transmitted signal power at time t; Let be the transmit antenna gain at time t; Indicates the receiver antenna gain; The distance between the transmitter and the mirror reflection point; The distance between the receiving end and the point of reflection on the mirror; The autocorrelation function of the pseudo-random code; The bistatic radar cross section per unit area is the normalized bistatic radar cross section (NBRCS). The normalized bistatic radar cross section The calculation formula is: in, For the calibrated signal power, , For atmospheric loss correction. For instrument loss correction, P represents the satellite transmitted signal power; This represents the transmitter antenna gain at the reflection point. This represents the receiver antenna gain at the reflection point. For distance loss from the transmitter to the sea surface and from the sea surface to the receiver, , For code correlation function and Doppler frequency shift function; The calculation of the normalized integral delay waveform leading edge slope (LES) is as follows: First, the integral delay waveform (IDW) is obtained using the delayed Doppler plot, and the calculation expression is: in, The number of Doppler frequencies is the unit of frequency. The correlation power of the sea surface scattered signal; The slope of the leading edge of the normalized integral delayed waveform, LES, is calculated using the slope equation. in, This represents the time delay value of the integral delay waveform IDW. The amplitude value of the integral delay waveform IDW. The number of sampling points used to calculate the slope of the front edge.

4. The method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R according to claim 1, characterized in that, The data sources for the dataset described in step 1 include: satellite observations and meteorological element information; The satellite observations include the longitude of the specular reflection point, the latitude of the specular reflection point, the incident angle of the specular reflection point, the distance from the transmitter to the specular reflection point, the distance from the receiver to the specular reflection point, the receiver gain at the specular reflection point, the delayed Doppler plot, the normalized bistatic radar cross section, and the slope of the leading edge of the delayed Doppler plot. The meteorological information includes wind speed and rainfall.

5. A method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R according to claim 4, characterized in that, Step 1 involves constructing the dataset, which includes: Satellite observations are processed, including removing missing values, removing data from areas other than water bodies, removing data with antenna gain below a threshold, and removing data with range correction gain below a threshold. The formula for calculating the range correction gain (RCG) is as follows: in, This represents the distance from the transmitter to the point of reflection on the mirror. This represents the distance from the receiver to the point of reflection on the mirror. This represents the receiver antenna gain at the point of specular reflection. Then, the reference wind speed and rainfall data from ERA5 are transformed in resolution using spatiotemporal bilinear interpolation to correspond with satellite observation data in time and space, forming a dataset for model training.

6. The method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R according to claim 1, characterized in that, In step 2: The fusion weight expression is: in, For the first The distance factor at time is used to fuse the weights. For the first The weighting of rainfall factors at different times. For the first The characteristics of attention at a distance in time For the first The characteristics of attention during rainfall at any given moment To fuse weights and generate the weight matrix of the network, To fuse the weights, the bias vector of the network is generated. For feature vector concatenation operations, For the first The hidden state output obtained by the LSTM network after processing the satellite observation feature parameters at time 1; softmax represents the activation function. The final expression for the wind speed inversion result is: in, For the first The final wind speed forecast at that moment. This is a linear transformation layer that maps the fused features to wind speed values.

7. A method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R according to claim 6, characterized in that, The distance attention head takes the height sequence as its main input and maps the geometric height information into a high-dimensional feature representation, namely the distance attention feature, through a distance encoder; the distance attention calculation expression is: in, This is the distance attention weight matrix. , , These are the query matrix, key matrix, and value matrix for distance attention, respectively. , , These are the weight matrices for generating the query matrix, the key matrix, and the value matrix, respectively, based on the distance to the attention head. For the first The hidden state output obtained by the LSTM network after processing the satellite observation feature parameters at a given time. For the first The distance at time step encodes the high-dimensional feature vector. For distance attention output, The dimension of the key vector; The rainfall attention head takes the rainfall intensity sequence as input and adaptively adjusts the weights of each frequency component based on the rainfall intensity to obtain the rainfall attention features: The rainfall attention calculation expression is: in, This is the rainfall attention weight matrix. , , These are the query matrix, key matrix, and value matrix for rainfall attention. , , These are the weight matrices for generating the query matrix, the key matrix, and the value matrix, respectively, for the rainfall attention head. For the first Rainfall at any given time is encoded into a high-dimensional feature vector. The output is the rainfall attention feature. is the dimension of the key vector.

8. A method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R according to claim 1, characterized in that, Step 3-1 is as follows: Step 3-1-1: The Kriging interpolation model for the boundary adaptive mesh is a targeted improvement on the classic Kriging interpolation model. Its basic expression is: in, For position interpolated predicted wind speed at the location, For the first Observation points The measured wind speed at the location, For the first Kriging weight coefficients for each observation point The total number of observations participating in the interpolation is denoted by , where is the Lagrange multiplier. Used to satisfy the unbiased estimation condition; The systematic expression for weight calculation in the classic Kriging interpolation model is: Expanded to: in, for 3D semi-mutation function matrix, This is a vector of weight coefficients. For the target vector, For position and The semivariogram values ​​between , ; Step 3-1-2: Based on the classic Kriging interpolation model, physical constraints on wind speed variation are introduced to obtain a semi-variogram model that conforms to the spatial variation law of wind speed, expressed as: in, Distance The semivariogram value at that location, For the nugget effect, For the partial sill value, For range parameters, The spatial distance between observation points These are the physical constraint weighting coefficients. It is a constraint function based on physical laws.

9. A method for monitoring and early warning of strong wind disasters based on BeiDou GNSS-R according to claim 1, characterized in that, Step 4 is as follows: Step 4-1: Obtain the time-series extrapolated wind speed through LSTM multi-step forward prediction, expressed as: in, For the first Wind speed forecast at any given time For the first LSTM hidden state at time step To predict the time step; Step 4-2: Based on the predicted wind speed, calculate the rate of change of wind speed, expressed as: in, For the first The rate of change of wind speed at any given time. For the first Wind speed acceleration at any given moment; Step 4-3: Record the wind speed values ​​output by the wind speed inversion model from Step 2-3 at consecutive time points. Use the sliding window coefficient of variation method to determine the stability of the wind speed inversion. The coefficient of variation within the time window is obtained based on the wind speed values ​​at multiple consecutive time points. The expression for calculating the sliding window coefficient of variation is: in, For the first Coefficient of variation at time 1 The standard deviation of wind speed within the sliding window. This represents the average wind speed within the sliding window. The length of the sliding window. For the first The actual wind speed inversion value at any given moment; Step 4-4: Combine the wind speed and wind speed change rate predicted by wind speed time series extrapolation with the coefficient of variation of the sliding time window to conduct a coordinated strong wind warning.

Images