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

Abstract

The invention discloses a Beidou GNSS-R-based strong wind disaster monitoring and early warning method, and the method comprises the steps: 1, enabling an unmanned plane to carry GNSS-R receiving equipment to be located at the periphery of a strong precipitation region, carrying out the wind speed detection of an internal water surface, and carrying out the data calculation of an obtained satellite observation quantity; step 2, constructing an LSTM wind speed inversion model fused with two-factor decoupling, performing inversion prediction on the water surface wind speed of a reflection point area by inputting the received satellite observed quantity and a resolving result and performing error correction, and obtaining future time step wind speed prediction based on a past time step wind speed inversion result; step 3, establishing a boundary adaptive grid and a dynamically updated Kriging interpolation model, obtaining wind speed data of a specific signal reflection point missing area, and assisting wind speed data supplement of the reflection point missing area; and 4, constructing a wind speed time sequence extrapolation and inversion stability coordinated double-effect strong wind early warning strategy, judging whether strong wind is about to appear or not, and performing early warning.

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 solve the above technical problems, the application discloses a strong wind disaster monitoring and early warning method based on Beidou GNSS-R, which comprises the following steps:

[0008] Step 1: the unmanned aerial vehicle carrying a GNSS-R receiving device is located at the periphery of a strong precipitation area to carry out wind speed detection operation on the water surface inside the area, data solving (GPS data solving is to obtain the position coordinates of a ground point by using the signal provided by a GPS satellite, receiving the signal by a receiver, and performing a series of complex mathematical processing) is performed on the obtained satellite observation data, and a data set is constructed;

[0009] Step 2: a LSTM wind speed inversion model fusing double-factor decoupling is constructed, the satellite observation data and the solving result are inputted and error correction is performed, the water surface wind speed of the reflection point area is inversely predicted, the wind speed at a future time step is predicted based on the wind speed at a past time step, for example, the wind speed at a future 15 minutes or 30 minutes is predicted according to the wind speed at a past 30 minutes;

[0010] Step 3: a boundary adaptive grid and a dynamic updating Kriging interpolation model are established, the water surface wind speed of the reflection point area is used for spatial adaptive interpolation, the expansion of the inversion area and the joint of the wind speed between the discrete reflection point areas are realized, the wind speed data of a specific signal reflection point missing area is obtained, and the wind speed data supplement of the reflection point missing area is assisted;

[0011] Step 4: a double-effect strong wind early warning strategy of wind speed time sequence extrapolation and inversion stability cooperation is constructed, the wind speed change trend is obtained by the future time step wind speed prediction provided by the LSTM wind speed inversion model, the variation coefficient of multiple continuous inversion wind speed values output by the LSTM wind speed inversion model is combined, and it is judged whether strong wind will occur and early warning is performed.

[0012] The data solving in step 1 is specifically as follows:

[0013] The GNSS-R receiver analyzes the characteristic change of the L-band signal emitted by the GNSS satellite after being reflected by the water surface, realizes the inversion of the water surface physical parameter, obtains the power distribution of the reflected signal, that is, the delay Doppler figure feature, and obtains the variable parameters related to the water surface roughness: the normalized bistatic radar scattering cross section NBRCS and the normalized integral delay waveform front slope LES, and the reflection point position information, signal quality information and other parameters are collectively used as the input of the LSTM wind speed inversion model.

[0014] The related power of the delay Doppler figure obtained by solving the original satellite signal is expressed as:

[0015]

[0016] Wherein, τ represents time delay; f krepresents the reflected signal power value of the selected delay Doppler frequency shift corresponding delay Doppler map (DDM) area; f represents the Doppler frequency shift; sinc(f) is the Doppler frequency shift function; T i represents the coherent integration time; λ represents the signal wavelength; P t represents the satellite transmitted signal power at t time; G t is the transmitted antenna gain at t time; G r represents the receiving antenna gain; R t is the distance from the transmitting end to the specular reflection point; R r is the distance from the receiving end to the specular reflection point; Λ(τ) is the pseudo-random code autocorrelation function; σ0 is the unit area bi-static radar scattering cross section, i.e. the normalized bi-static radar scattering cross section NBRCS, which reflects the sea surface scattering efficiency;

[0017] The normalized bi-static radar cross section σ0 calculation formula is:

[0018]

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

[0020] The calculation of the normalized integrated delay waveform front slope LES is as follows: the steepness of the delay Doppler map waveform front is described, and the distribution characteristics on the scattering surface represent the energy intensity distribution characteristics of the scattering signal. First, the integrated delay waveform IDW is obtained by using the delay Doppler map (DDM), and the calculation expression is:

[0021]

[0022] Wherein, N is the frequency unit Doppler number, and y is the sea surface scattering signal correlation power;

[0023] The normalized integrated delay waveform front slope LES is further calculated by the slope equation:

[0024]

[0025] Wherein, τ my is the time delay value of the integral delay waveform IDW m M is the number of sampling points participating in the calculation of the leading edge slope.

[0026] The data source of the data set in step 1 includes satellite observations and meteorological element information;

[0027] The satellite observations include mirror reflection point longitude, mirror reflection point latitude, mirror reflection point incident angle, transmitter to mirror reflection point distance, receiver to mirror reflection point distance, mirror reflection point receiver gain, delay Doppler map (DDM), normalized double base radar scattering cross section (NBRCS), and delay Doppler map leading edge slope (LES);

[0028] The meteorological element information includes wind speed and rainfall as reference values.

[0029] The construction of the data set in step 1 includes:

[0030] Processing the satellite observations includes removing missing values, removing data outside water bodies, removing data with antenna gain below a threshold value of 0, and removing data with range correction gain below a threshold value of 10, where the calculation formula of the range correction gain RCG is as follows:

[0031]

[0032] where, denotes the distance from the transmitter to the mirror reflection point, denotes the distance from the receiver to the mirror reflection point, denotes the receiver antenna gain at the mirror reflection point;

[0033] Then, the reference wind speed and rainfall data from ERA5 (ECMWF Reanalysis v5) are resolution-converted using spatio-temporal bilinear interpolation, corresponding to the satellite observation data in time and space, to form a data set for model training.

[0034] Step 2 specifically includes:

[0035] Step 2-1, using the LSTM algorithm to capture the time and spatial continuity of the reflection point trajectory, combining the fact that water surface wind speed also has continuity, training an LSTM basic wind speed inversion model, and obtaining the preliminary predicted wind speed inversion value by inputting the satellite observation feature parameters;

[0036] Step 2-2, calculating the double-factor decoupling error correction quantity for important influencing factors (distance and rainfall) in the wind speed inversion process, and obtaining the distance attention correction weight value and the rainfall attention correction weight value:

[0037] In the case of a certain satellite signal incidence angle, the signal reflection point information of different horizontal distances can be obtained by changing the flight height of the UAV, that is, the water surface wind speed condition of a farther distance can be detected by rising the height, but with the increase of the signal propagation path and the influence of rainfall, systematic deviation appears in the inversion algorithm.

[0038] Step 2-3, the fusion weight generation network takes the connection vector of the original LSTM feature, the distance attention feature and the rainfall attention feature as input, generates normalized fusion weight through nonlinear mapping, and outputs the final wind speed prediction result;

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

[0040] The fusion weight expression is:

[0041]

[0042] Among them, is the distance factor fusion weight at the t time, is the rainfall factor fusion weight at the t time, is the distance attention feature at the t time, is the rainfall attention feature at the t time, W f is the weight matrix of the fusion weight generation network, b f is the bias vector of the fusion weight generation network, [;;] is the feature vector connection operation, H t is the hidden state output obtained by the LSTM network processing satellite observation feature parameters at the t time;

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

[0044]

[0045] Among them, is the final wind speed prediction value at the t time, Linear is a linear transformation layer, which maps the fusion feature to a wind speed value.

[0046] By establishing the time sequence LSTM algorithm and adding double-factor error correction, the wind speed inversion precision is improved, which ensures that the UAV platform can obtain reliable wind speed information when located at different observation positions, and also provides technical support for future time step wind speed prediction.

[0047] The double-factor decoupling error correction in step 2-1 adopts a parallel double-head attention architecture, and a distance attention head and a rainfall attention head are respectively designed;

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

[0049]

[0050] wherein, is a distance attention weight matrix, are a query matrix, a key matrix and a value matrix of the distance attention respectively, are a weight matrix for generating the query matrix, a weight matrix for generating the key matrix and a weight matrix for generating the value matrix of the distance attention head respectively, H t is a hidden state output obtained by processing the satellite observation feature parameter through the LSTM network at the t th moment, is a distance encoding high-dimensional feature vector at the t th moment, is a distance attention output, d k is the dimension of the key vector;

[0051] The rainfall attention head takes the rainfall intensity sequence as the input, and adjusts the weight of each frequency component according to the rainfall intensity to obtain a rainfall attention feature: the rainfall attention calculation expression is:

[0052]

[0053] wherein, is a rainfall attention weight matrix, are a query matrix, a key matrix and a value matrix of the rainfall attention respectively, are a weight matrix for generating the query matrix, a weight matrix for generating the key matrix and a weight matrix for generating the value matrix of the rainfall attention head respectively, is a rainfall encoding high-dimensional feature vector at the t th moment, is a rainfall attention feature output, d k is the dimension of the key vector.

[0054] Because GNSS-R technology adopts passive receiving satellite signal mode, it cannot guarantee sufficient signal reflection points in a specific observation area at all times, so a space-time interpolation model needs to be established to supplement the wind speed data in the area where some signal reflection points are missing. Then, the unmanned aerial vehicle platform can be used to quickly detect different water areas, obtain as many regional water surface wind speed inversion results as possible by changing different observation positions, and assist in supplementing the wind speed data in the area where some signal reflection points are missing. The "missing" refers to the fact that there are no satellite signal reflection points in some areas at this moment, so the wind speed cannot be directly obtained through "signal inversion", for example, the object of a certain task is a lake, but at a certain moment, only the southwest corner of the lake has a reflection point, so the northeast corner of the lake is the "area where some signal reflection points are missing". The supplement is realized through the processes of steps 3-1 and 3-2, which can be simply understood as "data interpolation under certain rules" (insert new data into the area without data according to the characteristics of the existing data, and the new data does not depend on the reflection point inversion)

[0055] By establishing a boundary adaptive grid Kriging interpolation model, based on the classic Kriging interpolation theory, considering the irregular characteristics of the water body boundary, the physical characteristics of the wind field, and the confidence and high-frequency update characteristics facing the rescue demand, the boundary adaptive grid division mechanism and physical constraint conditions are introduced, and the uncertainty quantification and adaptive update plate are added, so as to obtain the optimal spatial coverage effect with the least observation cost.

[0056] Step 3 is specifically:

[0057] Step 3-1, on the basis of the classic Kriging interpolation model, the physical condition constraint of wind speed change is introduced to construct a semi-variogram function model that conforms to the spatial variation law of wind speed;

[0058] Step 3-2, construct a boundary adaptive grid system, the grid is the basic spatial unit of wind speed interpolation, a certain wind speed value in each grid represents the wind speed in this grid area, which is expressed as:

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

[0060] Where, G a is the adaptive grid system, T k is the kth triangular grid unit, P b is the set of water boundary constraint points, P o is the set of observation point positions, and ρ k is the grid density parameter of the kth area;

[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 grid confidence interval calculation expression is:

[0072]

[0073] CI α (s0) is a confidence interval with a confidence level of (1-α) at position s0, is a Kriging prediction value at position s0, z α / 2 is a standard normal distribution two-sided critical value, σ k (s0) is a Kriging prediction standard deviation at position s0;

[0074] Step 3-5, after obtaining the adaptive grid Kriging wind speed interpolation result and confidence, real-time weight updating and adjustment are performed, and the expression is:

[0075]

[0076] wherein, is the updated weight coefficient at the t+1 moment, is the weight coefficient at the t moment, η is a learning rate parameter, and the value range is 0.01 to 0.05, L t is a total loss function containing prediction error and physical constraint;

[0077] And further spatial correlation adaptive updating is performed to obtain a semivariogram function that changes synchronously with the new inversion wind speed data, so that the grid interpolation wind speed can be kept at a high frequency, and the expression is:

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

[0079] wherein, γ u (h) is the updated semivariogram function, γ o (h) is the semivariogram function before updating, γ e (h) is an empirical semivariogram function based on new observation data, and ω is an updating weight used to control the fusion proportion of new and old information, and the value range is 0 to 1.

[0080] The traditional Kriging interpolation algorithm can well capture the spatial correlation of meteorological elements. In view of the actual demand of unmanned aerial vehicle platform for inland water surface wind speed inversion, the interpolation grid is improved and physical constraint is added, or the actual wind speed data is filled in to realize wind speed information acquisition of multi-target area, and the interpolation result confidence evaluation and adaptive updating module is added to provide more real decision-making assistance for rescue personnel.

[0081] Step 3-1 is specifically:

[0082] Step 3-1-1, the Kriging interpolation model of the boundary adaptive grid is improved around the classic Kriging interpolation model, and its basic expression is:

[0083]

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

[0085] The weight solving system expression in the classic Kriging interpolation model is:

[0086] K·λ=k0

[0087] is expanded as:

[0088]

[0089] wherein, K is a (n+1)×(n+1) dimensional semi-variogram function matrix, λ is a weight coefficient vector, k0 is a target vector, γ(s i ,s0) is the semi-variogram value between positions s i and s0;

[0090] Step 3-1-2, on the basis of the classic Kriging interpolation model, the physical condition constraint of wind speed change is introduced, and a semi-variogram function model conforming to the spatial variation law of wind speed is obtained, and the expression is:

[0091]

[0092] wherein, γ(h) is the semi-variogram 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 weight coefficient, and Φ p (h) is a constraint function based on physical law.

[0093] Step 4 is specifically:

[0094] In the face of rapidly changing wind speed and possible strong wind conditions, under the condition of obtaining spatiotemporal continuous wind speed information, the future wind speed and change trend are perceived by forward prediction through past time step wind speed change, combined with the inherent characteristics of the accuracy of wind speed increase in the GNSS-R high wind speed inversion process, the sliding window is used to calculate the variation coefficient of the inversion result, and a double-effect strong wind warning strategy of wind speed time series extrapolation and inversion stability cooperation is constructed, and effective gale disaster warning is provided;

[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 is 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] DDM Delay-Doppler Map

[0110] NBRCS Normalized Bistatic Radar Cross Section

[0111] LES Leading Edge Slope

[0112] LSTM Long Short-Term Memory

[0113] IDW Integral Delay Waveform

[0114] RCG Range Corrected Gain

[0115] ERA5 ECMWF Reanalysis v5

[0116] Beneficial effects:

[0117] Technical level:

[0118] 1. Breakthrough traditional monitoring limitations: Through the innovative way of unmanned aerial vehicle carrying GNSS-R receiving equipment, the technical difficulties such as limited observation distance, visual obstruction and rainfall influence of traditional anemometer are overcome, and multi-direction, long-distance and all-weather water surface wind speed detection is realized.

[0119] 2. High-precision wind speed inversion technology: The LSTM wind speed inversion model is constructed by fusing double-factor decoupling, and through the double error correction of distance attention and rainfall attention, the wind speed inversion precision and stability under strong rainfall observation environment are significantly improved.

[0120] 3. Enhanced spatial coverage: The boundary adaptive grid and dynamic updating Kriging interpolation model are established, the wind speed data in the missing area of signal reflection point is supplemented, the utilization rate of satellite observation data is improved, the monitoring coverage is expanded, and the spatial continuity of wind speed information is provided.

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

[0122] Application level:

[0123] 1. Solve the demand for emergency rescue: In the case of rainstorm, the technical means of safe and rapid acquisition of meteorological information in the accident area is provided for the emergency situation of ship distress and personnel falling into the water, which provides key support for rescue decision-making.

[0124] 2. Improve monitoring flexibility: It has the ability to flexibly adjust the observation area, and can track and monitor the multi-target area and spatiotemporal continuity according to actual needs, and adapt to complex and changing rescue needs.

[0125] 3. Reduce monitoring cost and risk: Avoid the high cost of laying a large number of observation stations or increasing the number of guards, reduce the risk of damage to instruments and equipment in rainstorm disasters, and realize non-contact safe monitoring.

[0126] 4. Wide application prospect: It can be applied to the rescue of falling into the water, ship overturning, and the safety protection of various water activities such as river ferry and scenic sightseeing, and has good popularization and application value. BRIEF DESCRIPTION OF DRAWINGS

[0127] Figure 1 The algorithm flowchart of the present application. DETAILED DESCRIPTION

[0128] A strong wind disaster monitoring and early warning method based on Beidou GNSS-R, comprising:

[0129] Step 1: The unmanned aerial vehicle carrying GNSS-R receiving equipment is located outside the strong precipitation area to carry out wind speed detection operation of the water surface inside it, and the satellite observation data obtained is solved to construct a data set;

[0130] Step 2: Construct a LSTM wind speed inversion model that fuses double-factor decoupling, input the received satellite observation data and the solving result and correct the error, and predict the water surface wind speed of the reflection point area, and obtain the future time step wind speed prediction based on the past time step wind speed inversion result;

[0131] Step 3: Establish a boundary adaptive grid and a dynamically updated Kriging interpolation model, use the water surface wind speed of the reflection point area for spatial adaptive interpolation, realize the expansion of the inversion area and the joint of the wind speed between each discrete reflection point area, to obtain the wind speed data of the missing area of the specific signal reflection point, and assist the wind speed data supplement of the missing area of the reflection point;

[0132] Step 4, construct a double-effect strong wind warning strategy that coordinates wind speed time series extrapolation and inversion stability, obtain the wind speed change trend through the future time step wind speed prediction provided by the LSTM wind speed inversion model, and combine the variation coefficient of multiple continuous inversion wind speed values output by the LSTM wind speed inversion model to determine whether strong wind will occur and issue a warning.

[0133] The data solution of step 1 is specifically:

[0134] The GNSS-R receiver analyzes the characteristic changes of the L-band signal reflected by the water surface after being transmitted by the GNSS satellite to realize the inversion of the water surface physical parameters; based on the correlation between water surface roughness and wind speed, the power distribution of the reflected signal, i.e. the delay Doppler figure feature, is solved to obtain the variable parameters related to the water surface roughness: normalized bistatic radar cross section NBRCS and normalized integral delay waveform front slope LES, and then the reflection point position information and signal quality information.

[0135] The related power of the delay Doppler figure solved from the original satellite signal is represented as:

[0136]

[0137] Wherein, τ represents the time delay; f k represents the reflected signal power value of the selected time delay Doppler shift corresponding to the delay Doppler figure region; f represents the Doppler shift; sinc(f) is the Doppler shift function; T i represents the coherent integration time; λ represents the signal wavelength; P t represents the satellite transmission signal power at t; G t is the transmission antenna gain at t; G r represents the receiving antenna gain; R t is the distance from the transmission end to the mirror surface reflection point; R r is the distance from the receiving end to the mirror surface reflection point; Λ(τ) is the pseudo-random code autocorrelation function; σ0 is the unit area bistatic radar cross section, i.e. the normalized bistatic radar cross section NBRCS;

[0138] The calculation formula of the normalized bistatic radar cross section σ0 is:

[0139]

[0140] Wherein, p g,τ,f is the calibrated signal power, L a1 , L a2 is the atmospheric loss correction, l τ,f is the instrument loss correction, P represents the satellite transmission signal power, is the corresponding transmitter antenna gain at the reflection point, for the corresponding receiver antenna gain at the specular 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 the code correlation function and the Doppler shift function;

[0141] The calculation of the normalized integrated delay waveform front slope LES is as follows: first, the integrated delay waveform IDW is obtained by using the delay Doppler figure, and the calculation expression is:

[0142]

[0143] Wherein, N is the number of Doppler frequency units, y is the correlation power of the sea surface scattering signal;

[0144] Further, the normalized integrated delay waveform front slope LES is calculated by the slope equation:

[0145]

[0146] Wherein, τ m is the time delay value of the integrated delay waveform IDW, y m is the amplitude value of the integrated delay waveform IDW, and M is the number of sampling points participating in the calculation of the front slope.

[0147] The data source of the data set in step 1 includes 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 incidence 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 of the specular reflection point, the delay Doppler figure, the normalized bistatic radar scattering cross section and the front slope of the delay Doppler figure.

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

[0150] The construction of the data set in step 1 includes:

[0151] The satellite observations are processed, including removing missing values, removing data in areas other than water, removing data with antenna gain below a threshold, and removing data with range correction gain RCG below a threshold, wherein the calculation formula of the range correction gain RCG is as follows:

[0152]

[0153] Wherein, represents the distance from the transmitter to the specular reflection point, represents the distance from the receiver to the specular reflection point, represents the receiver antenna gain at the specular reflection point;

[0154] Then the reference wind speed and rainfall data from ERA5 are converted in resolution by using spatio-temporal bilinear interpolation, and are corresponded with satellite observation data in time and space to form a data set for model training.

[0155] Step 2 is specifically:

[0156] In step 2-1, the LSTM algorithm is used to capture the time and space continuity of the reflection point trajectory, and the water surface wind speed also has the characteristic of continuity. An LSTM basic wind speed inversion model is trained, and the satellite observation feature parameters are input to obtain the preliminary predicted wind speed inversion value.

[0157] In step 2-2, the double-factor decoupling error correction quantity is calculated for the important influencing factors in the wind speed inversion process, and the distance attention correction weight and the rainfall attention correction weight are obtained.

[0158] In step 2-3, the fusion weight generation network takes the connection vector of the original LSTM feature, the distance attention feature and the rainfall attention feature as input, generates normalized fusion weight through nonlinear mapping, and outputs the final wind speed prediction result.

[0159] The fusion weight expression is:

[0160]

[0161] Wherein, is the distance factor fusion weight at the t-th moment, is the rainfall factor fusion weight at the t-th moment, is the distance attention feature at the t-th moment, is the rainfall attention feature at the t-th moment, W f is the weight matrix of the fusion weight generation network, b f is the bias vector of the fusion weight generation network, [;;] is the feature vector connection operation, H t is the hidden state output obtained by the LSTM network processing satellite observation feature parameters at the t-th moment;

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

[0163]

[0164] Wherein, is the final wind speed prediction value at the t-th moment, Linear is a linear transformation layer that maps the fusion feature to a wind speed value.

[0165] The double-factor decoupling error correction in step 2-1 adopts a parallel double-head attention architecture, and a distance attention head and a rainfall attention head are designed respectively.

[0166] Distance attention head takes height sequence as the main input, and maps the geometric height information into high-dimensional feature representation, i.e., distance attention feature, through distance encoder; the distance attention calculation expression is as follows:

[0167]

[0168] wherein, is the distance attention weight matrix, are the query, key and value matrices of the distance attention, respectively, are the weight matrices of the distance attention head, respectively, t is the hidden state output of the LSTM network after processing the satellite observation feature parameters at the t-th moment, is the distance encoding high-dimensional feature vector at the t-th moment, is the distance attention output, k is the dimension of the key vector.

[0169] The rainfall attention head takes the rainfall intensity sequence as the input, and adaptively adjusts the weight of each frequency component according to the rainfall intensity to obtain the rainfall attention feature; the rainfall attention calculation expression is as follows:

[0170]

[0171] wherein, is the rainfall attention weight matrix, are the query, key and value matrices of the rainfall attention, respectively, are the weight matrices of the rainfall attention head, respectively, is the rainfall encoding high-dimensional feature vector at the t-th moment, is the output of the rainfall attention feature, k is the dimension of the key vector.

[0172] Step 3 is specifically:

[0173] Step 3-1: On the basis of the classical Kriging interpolation model, the physical condition constraint of wind speed change is introduced to construct a semi-variogram function model conforming to the spatial variation law of wind speed.

[0174] Step 3-2: An adaptive grid system is constructed, and the grid is taken as the basic spatial unit for wind speed interpolation. A certain wind speed value in each grid represents the wind speed in the grid area, which is 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, where γ(0) is the semi-variogram value at zero distance, and z

[0186] Step 3-4, after realizing the spatial interpolation of wind speed, the kriging model can provide the characteristics of the prediction variance of each interpolation point, and the confidence interval is calculated based on the prediction variance to calculate the confidence, when the prediction uncertainty at the grid point is large, it indicates that the area needs to supplement observation, and provides reference for unmanned aerial vehicle to replace new observation point;

[0187] The grid confidence interval calculation expression is:

[0188]

[0189] where CI α (s0) is the confidence interval with a confidence level of (1-α) at position s, is the kriging prediction value at position s0, z α / 2 is the two-sided critical value of the standard normal distribution, σ k (s0) is the kriging prediction standard deviation at position s0;

[0190] Step 3-5, after obtaining the adaptive grid kriging wind speed interpolation result and the confidence, real-time weight updating and adjustment are carried out, and the expression is:

[0191]

[0192] where, is the updated weight coefficient at the t+1 moment, is the weight coefficient at the t moment, η is the learning rate parameter, and the value range is 0.01 to 0.05, L t is the total loss function containing the prediction error and physical constraint;

[0193] And further spatial correlation adaptive update is carried out to obtain the semi-variogram function which changes synchronously with the new inversion wind speed data, so that the grid interpolation wind speed can be kept at high frequency update, and the expression is:

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

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

[0196] Step 3-1 specifically:

[0197] Step 3-1-1, the Kriging interpolation model of the boundary adaptive grid is improved around the classic Kriging interpolation model, and the basic expression is:

[0198]

[0199] wherein, is the wind speed interpolation prediction value at position s0, V(s i is the measured wind speed value at the i-th observation point s i , λ i is the Kriging weight coefficient of the i-th observation point, n is the total number of observation points participating in interpolation, is the Lagrange multiplier, and μ is used to satisfy the unbiased estimation condition;

[0200] The weight solving system expression in the classic Kriging interpolation model is:

[0201] K·λ=k0

[0202] is expanded as:

[0203]

[0204] wherein, K is a (n+1)×(n+1) dimensional semi-variogram matrix, λ is a weight coefficient vector, k0 is a target vector, γ(s i ,s0) is the semi-variogram value between positions s i and s0;

[0205] Step 3-1-2, on the basis of the classic Kriging interpolation model, the physical condition constraint of wind speed change is introduced, and a semi-variogram model conforming to the spatial variation law of wind speed is obtained, and the expression is:

[0206]

[0207] wherein, γ(h) is the semi-variogram value at a distance h, C0 is the block effect, C1 is the partial base value, a is the range parameter, h is the spatial distance between observation points, β is the physical constraint weight coefficient, and Φ p (h) is a constraint function based on physical law.

[0208] Step 4 is specifically:

[0209] Step 4-1, the time series extrapolation wind speed is obtained by multi-step forward prediction of LSTM, and the expression is:

[0210]

[0211] wherein, is the wind speed prediction value at the t+k time, h tis the LSTM hidden state at the t time, N is a prediction time step;

[0212] Step 4-2, based on the predicted wind speed, the wind speed change rate is calculated, and the expression is:

[0213]

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

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

[0216] Step 4-3, record the wind speed value obtained by the wind speed inversion model output in step 2-3 at continuous time, and use the method of calculating the sliding window coefficient of variation to judge the stability of the wind speed inversion, according to the wind speed value at multiple continuous time, the wind speed variation coefficient in the time window is obtained, if the variation coefficient is larger, it reflects that the wind speed is changing rapidly, otherwise it indicates that the wind speed value fluctuates smaller, that is, more stable, the sliding window variation coefficient calculation expression is:

[0217]

[0218] Where, CV t is the variation coefficient at the t time, σ t is the standard deviation of wind speed in the sliding window, μ t is the mean wind speed in the sliding window, W is the length of the sliding window, V i is the actual wind speed inversion value at the i time;

[0219] Step 4-4, combined with the wind speed time series extrapolation predicted wind speed and wind speed change rate and sliding time window variation coefficient calculation result, strong wind warning is carried out.

[0220] Example 1:

[0221] This embodiment takes the strong wind warning monitoring of the lake sightseeing tourist ship in the natural scenic area as an example to illustrate the implementation process of the present application:

[0222] Step 1, satellite observation data acquisition and calculation, get the characteristic input of wind speed inversion model:

[0223] A multi-rotor unmanned aerial vehicle equipped with GNSS-R receiving equipment is selected, the flight height is set to 200 meters, and it is located in the safety area outside the target lake area rainfall belt. The target monitoring area is an elliptical lake with an area of about 15 square kilometers, and there are 3 tourist sightseeing ships operating in the lake area.

[0224] GNSS-R receiver receives satellite direct and reflected signals in the area, obtains delay Doppler map through raw data solution, and extracts key parameters, including mirror reflection point coordinates, incident angle, normalized bistatic radar scattering cross section NBRCS value, normalized integral delay waveform front slope LES value, and distances of reflection point from signal transmitter and receiver respectively.

[0225] Step 2, fusion of double-factor decoupling LSTM wind speed inversion according to satellite observation data, to obtain the inversion wind speed at the signal reflection point:

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

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

[0228] Step 3, adaptive grid division of lake surface water area, and dynamic updating of improved Kriging model wind speed interpolation, to expand the lake surface wind speed monitoring range and improve the spatiotemporal resolution:

[0229] Based on the irregular boundary characteristics of the lake area, an adaptive grid system containing 458 triangular grid units is constructed. The basic grid density p0 is set to 0.025, the observation point distance attenuation coefficient a is set to 0.18, and the boundary complexity weight β c is set to 0.35.

[0230] Since the signal reflection point is missing in the northeast corner of the lake area at the current time, and the reflection is mainly concentrated in the middle of the lake area, the improved Kriging interpolation model is used to estimate and assign the grid wind speed in this area. The interpolation result shows that the wind speed in this area is 8.1 m / s, the confidence interval is [7.6, 8.6] m / s, and the confidence level is 92%.

[0231] Step 4, strong wind warning according to wind speed change:

[0232] The sliding window time length is set to 10 minutes, and the LSTM time series model is used to predict the wind speed change trend in the next 30 minutes. The prediction result shows that: 15 minutes later, the wind speed will reach 12.1 m / s, 30 minutes later, the wind speed will reach 15.7 m / s, and the wind speed change rate is 0.24 m / s 2.

[0233] At the same time, the wind speed variation coefficient CV in the current 10-minute window is 0.12, which exceeds the set threshold value 0.1, indicating that the wind speed is changing rapidly. Combining the time series extrapolation result and the variation coefficient analysis, the system predicts that strong wind weather (≥12 m / s) will occur within 30 minutes, and according to the geographical environment of the natural scenic spot and the wind and wave resistance level of the sightseeing ship, the warning level is formulated, including the first warning (wind speed ≥18 m / s, ship overturning probability is very large, need to immediately dock nearby or return to the docking port), the second warning (wind speed ≥15 m / s, ship overturning probability is larger, need to immediately organize the formation to return to the docking port in order), the third warning (wind speed ≥12 m / s, ship has the risk of overturning, limit the distance of ship from the shore), the unmanned aerial vehicle automatically sends the second warning signal to the scenic spot management duty room according to the wind speed inversion monitoring result.

[0234] Warning execution result:

[0235] The system sends warning information to the lake area management center and three tourist ships at the same time: "The wind speed in the lake area is expected to reach 15.7 m / s within 30 minutes, and the ship is required to immediately organize the formation to return to the port to avoid the wind." After receiving the warning, the ship returns to the port in time, successfully avoiding the safety accident that may be caused by strong wind weather.

[0236] The actual verification shows that the wind speed prediction error of this warning is less than 1.2 m / s, and the warning time is advanced by 28 minutes, which provides sufficient time guarantee for the safe evacuation of the ship.

[0237] The present application provides a strong wind disaster monitoring and warning method based on Beidou GNSS-R, and there are many methods and ways to realize this technical solution. The above description is only the preferred embodiment of the present application, and it should be pointed out that for ordinary technical personnel in this technical field, without departing from the principle of the present application, some improvements and refinements can be made, which should be regarded as the protection scope of the present application. The components not explicitly described in the embodiment can be realized by using existing technology.

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 an LSTM wind speed inversion model that integrates two factors and decouples them. By inputting the received satellite observations and the solution results and correcting the 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. 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. 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 obtained by solving the original satellite signal is expressed as follows: 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 R represents the receiver antenna gain; 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. The formula for calculating the normalized bistatic radar cross section σ0 is: 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; 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: Where N is the number of frequency units Doppler, and y is the correlation power of the sea surface scattered signal; The normalized integral delay waveform leading edge slope (LES) is further calculated using the slope equation: 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.

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, Step 2 is as follows: 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: 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. 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 fusion weight expression is: 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 is obtained by the LSTM network after processing the satellite observation feature parameters at time t, and softmax represents the activation function. The final expression for the wind speed inversion result is: 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.

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 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; 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 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; 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. 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.

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 specifically involves: 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: G a ={T k ∣T k =Delaunay(P b ∪P o ,ρ k )} 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 For the k-th region, Delaunay represents the construction of triangular mesh cells; 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: ρ(s)=ρ0·exp(-α·d o (s))·(1+β c ·C b (s)) 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; 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, 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; 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, Let be the Kriging prediction variance at position s0, and γ(0) be the semivariogram value at zero distance; 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. The expression for calculating grid confidence intervals is: 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; 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, 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; 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: c u (h)=(1-ω)·γ o (h)+ω·γ e (h) 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.

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 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, 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; The systematic expression for weight calculation in the classic Kriging interpolation model is: K·λ=k0 Expanded to: 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; 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: 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.

10. 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 specifically involves: Step 4-1: Obtain the time-series extrapolated wind speed through LSTM multi-step forward prediction, expressed as: 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. Step 4-2: Based on the predicted wind speed, calculate the rate of change of wind speed, expressed as: D 2 V t+k =ΔV t+k -ΔV t+k-1 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; 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: 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 mean wind speed within the sliding window. i The actual wind speed inversion value at time i; 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