Near-real-time water surface height inversion method

Abstract

The invention relates to the technical field of surveying and mapping science, in particular to a near-real-time water surface height inversion method, which comprises the following steps of: acquiring a shared frequency signal of each navigation system, setting different time windows according to different azimuth angles, and obtaining reflected signal sequence fragments of each navigation system at different azimuth angles; and calculating the reflected signal sequence fragments by using a spectral analysis method or a reverse modeling method to obtain the height of the reflecting surface, and performing inversion by using the height of the reflecting surface to obtain the height of the water surface. The problems that an existing model method is complex, a large amount of data is needed for modeling, results are prone to being influenced by noise, and efficiency and accuracy cannot be considered at the same time are solved.

Description

Technical Field

[0001] The present invention relates to the field of surveying and mapping science and technology, and in particular to a near real-time water surface height inversion method. Background Art

[0002] GNSS-IR (Global Navigation Satellite System Interferometric Reflectometry) is a technology that uses the multipath interference characteristics of GNSS signals to invert surface parameters. The frequency characteristics of the reflected signal are closely related to the height of the reflecting surface. Therefore, the water surface height can be inverted by analyzing the characteristic frequencies in the signal.

[0003] GNSS multipath interferometry inversion provides an innovative means for surface parameter monitoring by turning harm into benefit, that is, converting interference signals into effective information. It is particularly suitable for low-cost monitoring of large-scale, continuous dynamic environments and can be widely used in hydrological and ocean monitoring.

[0004] Patent publication number CN113625312A quantifies delayed sea state deviations by constructing model delays and observation delays. Based on the BM4 model, two parameters, reflection angle and incident angle, independent of sea conditions, are introduced, and the coefficients of each parameter are regressed to construct a sea state deviation parameter model. Based on the sea state deviation parameter model, delayed sea state deviations are predicted. This method has a complex model, requires a large amount of data for modeling, and is susceptible to noise interference, making it difficult to simultaneously achieve both efficiency and accuracy. Summary of the Invention

[0005] In response to the shortcomings of existing methods, the present invention shortens the required observation time window and improves the real-time performance of inversion without sacrificing inversion accuracy, and is widely applicable to the fields of hydrology and ocean monitoring.

[0006] The technical solution adopted by the present invention is: a near real-time water surface height inversion method comprises the following steps:

[0007] Step 1: Acquire the shared frequency signal of each navigation system, set different time windows according to different azimuth angles, and obtain the reflected signal sequence fragments of each navigation system at different azimuth angles;

[0008] As a preferred embodiment of the present invention, the navigation system includes GPS, Galileo and BDS-3.

[0009] As a preferred embodiment of the present invention, setting different time windows at different azimuth angles includes:

[0010] Determine azimuth and elevation ranges;

[0011] Set the length and step size of the time window;

[0012] Extract the dSNR of the same carrier frequency of each navigation system;

[0013] Convert the SNR into linear units and remove the trend term to obtain the dSNR of the reflected signal sequence fragments.

[0014] Step 2: Use spectrum analysis or inverse modeling to integrate and calculate the fragments of the reflected signal sequence to obtain the height of the reflecting surface, and use the height of the reflecting surface to invert and obtain the water surface height;

[0015] As a preferred embodiment of the present invention, the spectrum analysis method includes:

[0016] The dSNR values are averaged within the set elevation angle window to obtain a recombined dSNR sequence;

[0017] Adaptively denoise the reconstructed dSNR sequence;

[0018] LSP analysis is used to extract the characteristic frequencies of the recombined dSNR sequence and convert them into water surface height.

[0019] As a preferred embodiment of the present invention, the reverse modeling method includes:

[0020] Amplitude normalization is performed on the fragments of the reflected signal sequence;

[0021] As a preferred embodiment of the present invention, amplitude normalization uses Hilbert transform to identify the signal envelope and then performs normalization processing.

[0022] The amplitude-normalized reflected signal fragments are recombined to obtain the recombined dSNR sequence, which is:

[0023]

[0024] Where, λ is the carrier wavelength; k is the damping coefficient; s is the roughness parameter of the reflecting surface; e is the satellite elevation angle; h is the height of the reflecting surface; C1 and C2 are amplitudes;

[0025] Perform nonlinear fitting on the amplitude, reflective surface roughness parameters, and reflective surface height to obtain the optimal solution for the reflective surface height;

[0026] The water surface height is inverted using the optimal reflecting surface height.

[0027] As a preferred embodiment of the present invention, the nonlinear fitting adopts the nonlinear least squares method.

[0028] As a preferred embodiment of the present invention, it also includes quality control of the height of the reflecting surface.

[0029] As a preferred embodiment of the present invention, a near real-time water surface height inversion system includes: a memory for storing instructions executable by a processor; and a processor for executing the instructions to implement a near real-time water surface height inversion method.

[0030] As a preferred embodiment of the present invention, a computer readable medium stores computer program code, and the computer program code implements a near real-time water surface height inversion method when executed by a processor.

[0031] Beneficial effects of the present invention:

[0032] 1. This invention has discovered for the first time that recombining the shared frequencies of different navigation systems in the same time window can effectively suppress noise and amplify effective signals, thereby improving the accuracy of water surface height inversion;

[0033] 2. The present invention designs and restructures the dSNR model to improve the real-time performance and time resolution of the inversion, and adapts to the needs of dynamic monitoring;

[0034] 3. Normalization processing enhances noise resistance and improves fitting stability;

[0035] 4. The method of the present invention can be widely used in high-frequency water level monitoring scenarios in various water areas such as oceans, lakes, and reservoirs. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 It is a logic block diagram of the near real-time water surface height inversion method of the present invention;

[0037] Figure 2 It is a near real-time inversion sliding time window;

[0038] Figure 3 is the dSNR sequence combined elevation angle window;

[0039] Figure 4 is the amplitude normalization of the G24 satellite dSNR sequence;

[0040] Figure 5 is the time window SNR and dSNR sequence of AT01 station;

[0041] Figure 6 It is the LSP analysis of the time window combined dSNR sequence of AT01 station;

[0042] Figure 7 is the time window SNR and dSNR sequence of station SC02;

[0043] Figure 8 It is the LSP analysis of the time window combined dSNR sequence of station SC02;

[0044] Figure 9It is the time series water surface height inversion using the near real-time spectrum analysis method at AT01 station;

[0045] Figure 10 It is the time series water surface height inversion using the near real-time spectrum analysis method at SC02 station;

[0046] Figure 11 is the parameter fitting of the normalized dSNR sequence in the time window of station AT01;

[0047] Figure 12 is the parameter fitting of the normalized dSNR sequence in the time window of station SC02;

[0048] Figure 13 It is the time series water surface height inversion of AT01 station using the near real-time inverse modeling method;

[0049] Figure 14 It is the time series water surface height inversion of SC02 station using the near real-time inverse modeling method;

[0050] Figure 15 It is the comparison of the inversion RMSE and the number of inversion points between the classical method and the near real-time method;

[0051] Figure 16 It is the distribution statistics of the time intervals between adjacent inversion points at stations AT01 and SC02;

[0052] Figure 17 It is SC02 station DOY1367:40-8:10dSNR fragment LSP analysis;

[0053] Figure 18 This is the azimuth distribution of DOY201 signals in different time windows at AT01 station;

[0054] Figure 19 It is the difference distribution of the inversion results of near real-time spectrum analysis and inverse modeling method in the same time window;

[0055] Figure 20 It is the dSNR sequence incomplete window spectrum analysis method and the inverse modeling method for water surface height inversion;

[0056] Figure 21 It is the inversion statistics of the data missing window of the spectrum analysis and inverse modeling method for stations AT01 and SC02. DETAILED DESCRIPTION

[0057] The present invention will be further described below in conjunction with the accompanying drawings and embodiments. This figure is a simplified schematic diagram, which only illustrates the basic structure of the present invention in a schematic manner, and therefore only shows the components related to the present invention.

[0058] Traditional spectrum analysis methods such as Lomb-Scraggle (LSP) rely on complete low-altitude angle observation data, which often takes 40 minutes to 1 hour and is subject to significant noise interference. Post-processing is required to remove outliers, affecting real-time and timeliness. The present invention discovers for the first time that recombining the shared frequencies of different navigation systems in the same time window can effectively suppress noise and amplify effective signals, thereby improving the accuracy of water surface height inversion.

[0059] like Figure 1 As shown, a near real-time water surface height inversion method includes the following steps:

[0060] Step 1: Acquire different navigation shared frequency signals and perform preprocessing;

[0061] In this embodiment, two stations, AT01 and SC02, which can simultaneously receive GPS, Galileo, and BDS-3 shared frequency signals, are selected. Near-real-time water surface height inversion is performed using the SNR data of the two shared frequencies of GPS, Galileo, and BDS-3 at AT01 and SC02. The signal characteristics are shown in Table 1. The signals from GPS L1, Galileo E1, and BDS-3 B1C all have a frequency of 1575.42 MHz and are combined to form the C1 signal group. The signals from GPS L5, Galileo E5a, and BDS B2a all have a frequency of 1176.45 MHz and are combined to form the C5 signal group.

[0062] Table 1 GPS, Galileo, and BDS-3 shared frequency signals and characteristics

[0063]

[0064] A sliding time window is used to intercept the dSNR fragments of the reflected signals from different systems and satellites, and then reassemble them to form a dSNR sequence within the window.

[0065] Different from the existing methods that mainly use complete low-altitude angle SNR data, the present invention first discovers and proposes to focus on the SNR data fragments of different satellites within a sliding time window and try to extract the reflection surface height from these fragments.

[0066] like Figure 2 As shown in the figure, taking GPS L1, Galileo E1 and BDS-3B1C as examples, the steps for extracting SNR data fragments in the sliding time window include:

[0067] 1. Determine the available azimuth and elevation angle ranges, i.e., filter the reflected signals from the water surface;

[0068] 2. Determine the time sliding window length and sliding step size;

[0069] 3. Extract SNR data of the same carrier frequency from the RINEX file;

[0070] 4. Select satellite tracks that meet the azimuth, elevation, and time window conditions; be careful to remove satellite tracks that are too short (e.g., less than 10 epochs);

[0071] 5. Convert the raw SNR data of each satellite into linear units and use a low-order polynomial to remove the trend term to obtain the dSNR sequence fragments of the reflected signals of different satellites in the time window.

[0072] The low-order polynomial has an order less than 3.

[0073] When extracting SNR data fragments, the selection of the time window length is crucial. Using a smaller time window will result in insufficient SNR data, leading to unreliable inversion results. However, using a larger time window is not conducive to improving the real-time performance of the inversion. The length of the time window needs to take into account the available water surface azimuth range of the station and the number of visible GPS, Galileo, and BDS-3 satellites, of which the available water surface azimuth range plays a decisive role. A larger range of available azimuth angles means more satellite tracks within each time window, allowing for a shorter time window. Generally speaking, for a station that can receive GPS, Galileo, and BDS-3 signals, when the available azimuth angle range is greater than 200°, a time window of approximately 20 minutes can be selected. When the available azimuth angle range is between 150° and 200°, a time window of approximately 30 minutes is required. For stations with an available azimuth angle range less than 150°, a time window of more than 40 minutes is required. The sliding step size of the time window can be selected based on the time resolution requirements of different application scenarios and is not limited.

[0074] Step 2: Use spectrum analysis or inverse modeling to integrate and calculate the fragments of the reflected signal sequence to obtain the height of the reflecting surface, and use the height of the reflecting surface to invert and obtain the water surface height;

[0075] The construction of the spectrum analysis method includes:

[0076] Existing methods process SNR sequence fragments of several time lengths separately, i.e., perform LSP analysis on each SNR sequence fragment. The present invention finds that in LSP analysis, the accurate identification of characteristic frequencies depends on the length of the signal. When the data is truncated, the inversion outliers will increase significantly, resulting in a decrease in inversion accuracy. The present invention also finds that when the signal is confined to a smaller time window, the reflected signals from the same frequency carrier have similar characteristics.

[0077] The reflected signal dSNR is usually expressed as:

[0078]

[0079] Where A is the reflected signal amplitude, is the phase, λ is the carrier wavelength, h is the height of the reflecting surface, and sine is the time variable that changes with the altitude angle.

[0080] The signal amplitude A part depends on the satellite transmission power and varies with different satellites; the phase It is a frequency-related quantity, which mainly depends on the electromagnetic characteristics of the reflecting surface. In each time window, the electromagnetic characteristics of the water surface are relatively fixed, and the phase of the dSNR sequence of different satellites can be considered to be Ignoring the change of water surface in the time window, h in the expression of dSNR signals of different satellites is also the same; this means that the dSNR signals with the same frequency carrier also have the same signal frequency in the time window; assuming that there are N segments of dSNR fragments with the same frequency carrier in a time window, it can be seen that the superposition of these N segments of signals can be understood as frequency Same, phase The same signal with different amplitude A is superimposed; therefore, the result is still a cosine signal with the same frequency; the amplitude of the recombined signal oscillates with time, but the frequency does not change. Using LSP analysis, the characteristic frequency can still be extracted, and then the water surface height can be inverted.

[0081] When using the spectrum analysis method, after obtaining the dSNR fragments of each satellite in the time window, the dSNR values are averaged in a certain altitude angle window. Then, based on the preliminary fitting results of the signal, some obvious outliers are eliminated to form the dSNR reconstructed signal in the time window, such as Figure 3 As shown in the figure, lines of different colors represent dSNR sequences from different satellites, and black scattered dots represent dSNR data points within the elevation angle window. This process smoothes the entire dSNR sequence, reduces the impact of edge noise, and highlights the main trend of the signal. However, it should be noted that an elevation angle window that is too wide will change the signal characteristics. For observations with a sampling rate of 15s, it is recommended to use an elevation angle window of 0.05°.

[0082] After obtaining the recombined dSNR sequence within the elevation angle window, an adaptive denoising method is used to separate the noise signal and weaken the influence of low-frequency and high-frequency noise on the characteristic frequency identification.

[0083] Then, LSP analysis is used to extract the characteristic frequency of the recombined dSNR within the time window, which is then converted into water surface height;

[0084] Among them, adaptive denoising uses wavelet decomposition.

[0085] Finally, certain quality control measures are adopted to eliminate outliers.

[0086] The selection of quality control measures refers to the classic GNSS-IR spectrum analysis method. The ratio of the main frequency power to the average background noise power P / N in the LSP and the spectrum amplitude peak A are used as quality control conditions. The time window that does not meet the conditions is regarded as an abnormal window and no inversion value is output.

[0087] Use the same process again to slide and process the next time window.

[0088] Build a near real-time inverse modeling approach, including:

[0089] As for the inverse modeling method, the study also considers the use of a time window to model the dSNR within the window. The difference is that instead of considering only the dSNR sequence of a single satellite as the modeling object, the dSNR fragments with the same frequency carrier within the time window are modeled as a whole to reconstruct the complete oscillation signal, thereby extracting enough oscillation waveforms to fit the optimal parameters. The inverse modeling method relies on the accurate fitting of the dSNR waveform. In order to better represent the oscillation of the signal, the signal model used is slightly different from that used in the spectrum analysis method.

[0090] The dSNR sequence oscillation will gradually decrease as the elevation angle increases. Therefore, a damping coefficient k is introduced, and the signal model is shown in formula (2):

[0091]

[0092] Where k = 2π / λ, which is the wave number of the GNSS signal; s is the roughness parameter of the reflecting surface, which is related to the amplitude of the water surface fluctuation.

[0093] Use C1 and C2 to replace the amplitude and phase in the original expression to ensure the stability of the inversion value;

[0094] The conversion relationship between C1 and C2 and the original phase and amplitude is:

[0095]

[0096] When using nonlinear least squares fitting, Strandberg et al. chose to share the C1 and C2 parameters for signals with the same carrier frequency over a long period of time (usually several days), and share the s parameter among signals of all frequencies. They then solved h for each satellite trajectory to obtain the water surface height at different time points. Therefore, it is essentially still a post-processing solution method. In order to obtain near real-time inversion results, the present invention adopts a different solution approach.

[0097] Within a time window, the dSNR signals from the GPS / Galileo / BDS-3 shared frequency carrier have the same signal frequency and phase. However, due to the different amplitudes of the reflected signals from different satellites, the amplitude of the combined dSNR signal will oscillate irregularly over time, which is not conducive to high-precision fitting of signal parameters. Here, we first perform amplitude normalization on the dSNRs of different satellites.

[0098] Hilbert Transform (HT) is a mathematical transformation in signal processing, which is widely used to analyze the envelope, phase and other information of the signal and effectively process non-stationary signals.

[0099] The dSNR fragment signal envelope of each satellite is identified using Hilbert transform, and the maximum amplitude is normalized;

[0100] Assume the original dSNR signal is x(t), and the Hilbert transform is Combining the original signal with its Hilbert transform yields the analytical signal z(t):

[0101]

[0102] Where j is the imaginary unit;

[0103] The amplitude of the analytical signal z(t) is the instantaneous amplitude A(t) of the signal. The instantaneous amplitude can be obtained from the modulus of the analytical signal:

[0104]

[0105] In order to reduce the noise interference at the edge of the signal, the average value of the first 1% of the maximum instantaneous amplitude of each trajectory is taken as the maximum amplitude for normalization; Figure 4 a represents the original dSNR (black solid line) and the instantaneous amplitude of the signal identified by Hilbert transform (red solid line); Figure 4 b is the dSNR sequence that has been amplitude normalized.

[0106] Within a time window, the normalized dSNR fragments have the same frequency, amplitude, and phase, that is, they have the same C1, C2, and h parameters. The reflector roughness parameter s is independent of the satellite and can be considered a fixed value within the time window. The normalized dSNR fragments are recombined and sorted according to the time variable sine to form a complete combined signal.

[0107] Then, the four unknown parameters C1, C2, h, and s are fitted by the nonlinear least squares method to obtain the optimal parameter solution, and the height h of the reflecting surface can be obtained, thereby obtaining the water surface height within the time window.

[0108] Inverse modeling is inherently an iterative approach, so the selection of initial parameters is crucial for nonlinear least-squares fitting. While iteration is insensitive to the C1, C2, and s parameters, it is particularly sensitive to the initial value of h. Inappropriate initial parameters can cause the parameters to converge to a local optimum, leading to inversion anomalies. Given the continuity of water surface changes within each time window, the inverted water surface height prior to the current time window can be used as the initial parameter.

[0109] Similar to the near-real-time spectrum analysis method, the near-real-time inverse modeling method also uses dSNR sequence fragments within the time window for water surface height inversion. The signal fragments after normalization share solution parameters, which provides favorable conditions for the realization of high-precision fitting.

[0110] Analysis of near-real-time water surface height inversion test results:

[0111] Using two sets of shared frequency signals of GPS, Galileo, and BDS-3: GPS L1, Galileo E1, BDS-3B1C and GPS L5, Galileo E5a, BDS-3B2a, near real-time inversion spectrum analysis method and inverse modeling method were tested at AT01 and SC02 stations, and compared with the classical inversion method.

[0112] Analysis of time window inversion results using spectrum analysis method:

[0113] The AT01 station's water surface azimuth range is 0°-220°. SNR data with elevation angles of 5°-13° are selected and used in a 20-minute time window to perform surface height inversion using spectrum analysis. During the inversion, a 0.05° elevation window is used, and the dSNR fragments within the window are averaged to smooth the entire dSNR sequence, highlighting the main trends of the signal and reducing the interference of edge noise. Figure 5 The SNR and dSNR sequences of GPS S1C, Galileo S1C, and BDS S1P signals in four different time windows at AT01 station on DOY 196, 2023 are shown; lines of different colors represent data from different satellites; on the left Figure 5 a, c, e, and g are the SNR data within the time window; Figure 5 b, d, e, and f are the dSNR sequences obtained after SNR detrending; Figure 5As an example, there are four satellites in the time window: G25, E32, C25, and C33. It can be seen that G25, E32, and C25 all have ascending elevation trajectory, and the SNR oscillation gradually decreases; C33 has a descending elevation trajectory, and the SNR oscillation gradually increases. The G25 S1C signal appears from the 397th epoch to the 440th epoch; the E32 S1C signal appears from the 379th epoch to the 440th epoch; the C25 S1P signal appears from the 361st epoch to the 440th epoch; the C33 S1P signal appears from the 420th epoch to the 440th epoch. The dSNR series obtained after detrending the SNR data of each signal are displayed in the right subplot according to the time variable sine. It can be seen that the amplitudes of the dSNR series of each satellite are slightly different, but the oscillation trends have shown a high degree of consistency. The black scattered points in the figure represent the window dSNR series formed by the combination.

[0114] right Figure 5 The four time windows in the dSNR sequence are combined for LSP analysis, and the results are as follows Figure 6 As shown in the figure, it can be seen that the LSP analysis results of the four windows all show a relatively obvious single peak, which can better extract the characteristic frequency of the reflected signal. The converted reflection surface heights are: 12.710m, 12.690m, 12.684m, 12.711m, and the corresponding water surface height inversion errors are: -0.255m, -0.175m, -0.118m, -0.149m, respectively.

[0115] like Figure 7 Compared with AT01, the water surface azimuth angle range of SC02 is smaller, 50°-240°; therefore, the time window used is slightly larger than that of AT01, with a window length of 30 minutes; the elevation angle range is limited to 5°-13°; the spectrum analysis method and data processing strategy are the same as those of AT01. Figure 7 The SNR and dSNR sequences of GPS S1C, Galileo S1C, and BDS S1P signals in four different time windows at SC02 station on DOY 136, 2024 are shown. Lines of different colors represent data from different satellites. The left sub-figures a, c, e, and g are the SNR data in the time window; the right sub-figures b, d, e, and f are the dSNR sequences in the time window; each time window contains rising and falling trajectories from different azimuths; the black scattered points in the figure represent the window dSNR sequence formed by the combination; LSP analysis is performed on the combined sequence, and the results are shown in Figure 8In the four time windows, the LSP analysis results show that the reflector heights are 6.039m, 4.845m, 4.402m, and 4.888m respectively; the inversion errors are -0.037m, +0.082m, +0.005m, and -0.083m respectively; the spectrum analysis method has achieved good inversion results in the test windows of AT01 and SC02 stations; it is also noted that in the time window of 7:40-8:10, the combined dSNR sequence of SC02 station is not continuous, such as Figure 8 Get f.

[0116] Analysis of time series inversion results: In order to further verify the effectiveness of the near-spectrum analysis method, the observation data of AT01 station from DOY196 to DOY215 in 2023, a total of 20 days, were used to perform time series water surface height inversion using the near-real-time spectrum analysis method for the C1 signal group formed by GPS S1C, Galileo S1C and BDS S1P, and the C5 signal group formed by GPS S5Q, Galileo S5Q and BDS S5Q; the inversion results of the C1 and C5 signal groups in each time window were averaged to form the combined inversion result of C1+C5; the time window length remained unchanged, with a time window of 20 minutes and a sliding step of 5 minutes; within the restricted elevation angle range, time windows with a combined dSNR sequence length less than 4° were regarded as invalid windows, and no inversion values were output; in addition, time windows with LSP analysis P / N>2.0 and A>2 were regarded as valid windows, and those that did not meet the conditions were regarded as invalid windows, and no inversion values were output.

[0117] The SC02 station used 15 days of observation data from DOY134 to DOY148 in 2024 to conduct a water surface height inversion test using the near-real-time spectrum analysis method for the C1, C5 and C1+C5 signal groups; a sliding window with a time window length of 30 minutes and a sliding step of 5 minutes was used; the reflecting surface height of the SC02 station was about 5.4 meters, which was lower than that of the AT01 station (12.6 meters). At the same time, there was less interference oscillation; therefore, a longer dSNR sequence length was required for the selection of the effective time window, and the time window with a combined dSNR sequence elevation angle length of less than 5° was regarded as an invalid window; in addition, the time window with LSP analysis P / N < 2.0 and A < 2 was regarded as an invalid window, and no inversion value was output for further quality control.

[0118] Figure 9 and Figure 10 The water surface height time series inversion results of the spectrum analysis method for stations AT01 and SC02 are presented respectively; Figure 9-10a is the comparison between the inversion results of the near-real-time spectrum analysis method for each signal group and the tide level observations; blue scatter points represent the inversion results of the C1 signal group, green scatter points represent the inversion results of the C5 signal group, and red scatter points represent the inversion results of the C1+C5 signal group. The tide level observations are represented by the black solid line; Figure 9-10 b is the number of time windows in which each signal group can output inversion values every day. It can be seen that using the near-real-time spectrum analysis method, the inversion results of each signal group at AT01 and SC02 stations are highly consistent with the tide level observations. Since the C5 signal group currently has less observation data than the C1 signal group, the number of time windows in which the C5 signal group can output inversion values every day is less than that of the C1 signal group. In almost all time windows in which the C5 signal group can provide valid inversion values, the C1 signal group can also output valid inversion results. Therefore, the number of valid time windows of the C1+C5 combination in the figure is very close to that of the C1 signal group.

[0119] Table 2 shows the number of valid windows, window efficiency, inversion root mean square error (RMSE), mean absolute error (Mea), and correlation coefficient R between inversion values and tidal observation values for the near-real-time spectrum analysis methods C1 signal group, C5 signal group, and C1+C5 signal group at AT01 and SC02 stations, which can normally output inversion values. The window efficiency in the table refers to the ratio of the number of valid time windows to the total theoretical number of windows. The time window of AT01 station is set to 20 minutes, the sliding step is 5 minutes, and there are 5700 theoretical time windows in 20 days. The time window of SC02 station is set to 30 minutes, the sliding step is 10 minutes, and there are 2130 theoretical time windows in 15 days.

[0120] It can be seen that at the AT01 station, the C1 signal group output a total of 5398 time windows that can output inversion results during the experiment, with a window efficiency of 94.7%, and inversion RMSE, Mea and R of 0.147m, 0.122m and 0.953 respectively; the window success rate and inversion accuracy of the C5 signal group are slightly lower than those of the C1 signal group; the combined inversion of C1 and C5 signals has 5410 windows that can output inversion values, with a window efficiency of 94.9%, and inversion RMSE, Mea and R of 0.149m, 0.127m and 0.953 respectively. The available azimuth inversion at SC02 station is smaller than that at AT01 station, and the available SNR data in each time window is also less. The window efficiency of the near-real-time spectrum analysis method is slightly lower than that of AT01 station; but it can be seen that the window efficiency of the spectrum analysis method for the C1 signal group of SC02 station can reach 90.9%; the window success rate of the C5 signal group is slightly lower than that of the C1 signal group, which is 86.1%; in addition, the correlation coefficient R between the inversion results of each signal group using the spectrum analysis method at SC02 station and the tide level observation value is better than 0.99; comparing the inversion results of the two stations, it can be seen that the overall inversion accuracy of SC02 station is better than that of AT01 station; this may be because the tide gauge station corresponding to AT01 station is far away from the GNSS receiver (74km), while SC02 is closer to the tide gauge station (359m), so more accurate tide level observation values can be used as reference values.

[0121] Table 2 Statistics of time series water surface height inversion results using near real-time spectrum analysis method

[0122]

[0123] Analysis of time window inversion results of near real-time inverse modeling method:

[0124] For water surface height inversion using the near-real-time inverse modeling method, the average of the top 1% maximum amplitudes after Hilbert transform is taken as the maximum amplitude for normalization. When performing parameter fitting using the nonlinear least squares method, the initial values of C1 and C2 are set to 1, and the initial value of s is set to 0. Considering the continuity of water surface changes, the average of the inversion results of the three time windows before the current time window is used as the initial estimate of h to further control the fluctuation range of h and improve the fitting accuracy.

[0125] Figure 11 and Figure 12These are the inversion results of the four time windows of AT01 and SC02 using the inverse modeling method, and the time windows remain unchanged; the scattered points of different colors in the figure represent the dSNR sequences from different satellites; it should be noted that the dSNR sequences in the figure have been amplitude normalized; it can be seen that as the altitude angle decreases, the amplitude of each signal gradually decreases; the black solid line in the figure represents the fitted reflection signal, which can better reflect the oscillation trend of each signal fragment; from the parameter fitting results, it can be seen that the reflection surface heights of the four windows of AT01 station are 12.694m, 1 The heights of the reflecting surfaces of the four windows of SC02 station are 2.716m, 12.706m, and 12.689m, and the inversion errors are -0.239m, -0.201m, -0.140m, and -0.127m respectively; the heights of the reflecting surfaces of the four windows of SC02 station are 6.021m, 4.771m, 4.370m, and 4.775m, and the inversion errors are -0.019m, +0.156m, +0.037m, and +0.030m respectively; the water surface height inversion values of each window of AT01 station and SC02 station using the near real-time inverse modeling method are consistent with the spectrum analysis method.

[0126] Analysis of time series inversion results: The observation data of stations AT01 and SC02 remain unchanged, and time series water surface height inversion using the near-real-time inverse modeling method is performed for the C1 signal group and the C5 signal group. The inversion results of the C1 and C5 signal groups in each time window are averaged to form the combined inversion result of C1+C5. The time window lengths of the two stations remain unchanged. In terms of quality control, the length of the normalized dSNR sequence within the window is used as the selection criterion. At stations AT01 and SC02, windows with window combination dSNR sequence lengths greater than 70 and 60, respectively, are used as valid time windows.

[0127] Figure 13 and Figure 14 The water surface height time series inversion results of the near-real-time inverse modeling method for stations AT01 and SC02 are shown respectively. It can be seen that using the near-real-time inverse modeling method, the inversion results of each signal group also have a high degree of consistency with the tide level observations. Similar to the near-real-time spectrum analysis method, the number of time windows in which the C1 signal group of the two stations can output inversion values every day is higher than that of the C5 signal group.

[0128] Table 3 summarizes the inversion results of the near-real-time inverse modeling method for stations AT01 and SC02. It can be seen that when using the near-real-time inverse modeling method, the window efficiency of the C1 signal group at both AT01 and SC02 stations is greater than 90%. Due to the small amount of available SNR data, the window efficiency of the C5 signal group is slightly lower than that of the C1 signal group, and the inversion accuracy is also slightly lower than that of the C1 signal group. Similar to the results of the near-real-time spectrum analysis method, the correlation coefficients between the inversion values of each signal group at AT01 and the tide level observation values are all greater than 0.95, and the correlation coefficients of each signal group at SC02 are all greater than 0.99. Comparing Tables 2 and 3, it can be noted that at AT01 station with a larger water surface azimuth angle, the window success rate of the spectrum analysis method for the three signal groups is higher than that of the inverse modeling method. At SC02 station with a smaller water surface azimuth angle coverage, the inverse modeling method can output inversion results in more time windows than the spectrum analysis method.

[0129] Table 3 Statistics of time series water surface height inversion results using the near real-time inverse modeling method

[0130]

[0131] Comparative analysis between near-real-time inversion method and classical inversion method:

[0132] The water surface height inversion of the C1, C5, and C1+C5 signal groups at stations AT01 and SC02 was performed using the classic GNSS-IR water surface hyperspectral analysis method. This method uses complete low-elevation angle trajectories without any time window, and the middle moment of each trajectory is used as the inversion value corresponding to the time. The inversion results of the C1 signal group are obtained using the low-elevation angle trajectories of GPS L1, Galileo E1, and BDS3 B1C. The inversion results of the C5 signal group are obtained using GPS S5Q, Galileo S5Q, and BDS3 S5Q signals. The inversion results of the C1+C5 signal group are obtained using the six signals included in the C1 and C5 signal groups. During quality control, the inversion values with LSP analysis P / N>3 and A>5 were selected as valid trajectories, and no other post-processing quality control measures were taken.

[0133] Figure 15The inversion RMSE and number of inversion points of the classical method and the near-real-time method during the AT01 and SC02 experiments are shown. In terms of inversion accuracy, compared with the classical method, the inversion RMSE of each signal group of the near-real-time spectrum analysis method and the inverse modeling method at the AT01 and SC02 stations are lower than that of the classical method, and the near-real-time inversion method has better accuracy than the classical inversion method. Gholamrezaee et al. used the classical method to evaluate the inversion results of multiple signals at the AT01 station with a time resolution of 1 hour, and the inversion accuracy was about 0.24m-0.36m. The near-real-time inversion spectrum analysis method and the inverse modeling method at the AT01 station were significantly better than this accuracy, which may be due to the stricter quality control adopted in this embodiment. The results of the study were published in the journal Nature Communications. Observations from multiple frequencies of GPS, GLONASS, Galileo, and BDS at SC02 were post-processed (outliers were removed within a 1-hour window) to invert water height. Adaptive denoising was used to achieve an RMSE of 0.110 m for the multi-system, multi-frequency inversion. Both the near-real-time spectrum analysis method and the inverse modeling method for the C1+C5 combination at SC02 outperformed this accuracy. In terms of the number of inversion points, the near-real-time spectrum analysis method increased the number of inversion points for the C1, C5, and C1+C5 signal groups at AT01 by an average of 142.3% compared to the classical method. At SC02, the average increase in the number of inversion points for each signal group was 149.2%. The near-real-time inverse modeling method increased the number of inversion points for each signal group at AT01 by an average of 141.9%, while the increase for SC02 was 152.1%.

[0134] It is worth noting that the near-real-time inversion method outputs the inversion value in each time window, so the time interval between adjacent inversion points is fixed and depends on the sliding window step size; however, the classical method uses a single-trace inversion, so the inversion time interval between adjacent inversion points is not consistent; Figure 16 The distribution of time intervals between adjacent inversion points of the near-real-time spectrum analysis method, inverse modeling method and classical inversion method at AT01 and SC02 stations was statistically analyzed. It can be seen that more than 90% of the inversion time intervals at AT01 and SC02 stations are fixed at 5 minutes and 10 minutes, respectively. However, when the inversion method is used, at AT01 station, 66.8% of the inversion time intervals are greater than 5 minutes, and the maximum time interval between adjacent inversion points is 74 minutes. At SC02 station, 59.5% of the inversion time intervals are greater than 10 minutes, and the maximum time interval between adjacent inversion points is 146 minutes. From the above analysis, it can be seen that compared with the classical GNSS-IR method, the near-real-time spectrum analysis and inverse modeling methods can provide water surface height inversion results with better accuracy, stronger real-time performance and more uniform time resolution.

[0135] There are many reasons why the near-real-time inversion method can achieve high-precision inversion values. First, the near-real-time inversion method uses dSNR sequence fragments from different satellites in the time window. Each trajectory carries the same characteristic signal frequency, but is subject to different noise interference. When different trajectory fragments are combined, the influence of noise is dispersed, and the characteristic frequency signals are combined and amplified, which is more conducive to the identification of characteristic frequencies. Therefore, the near-real-time method has better noise resistance than the classical method. Figure 17 The LSP analysis results of the dSNR fragments of each satellite for the C1 signal group of station SC02 in the time window of DOY1367:40-8:10 are shown; different colored lines represent the LSP analysis results of different satellites; it can be seen that the reflection surface height extracted from each satellite fragment is quite different, so it is impossible to accurately determine the water surface height; the LSP analysis results of the recombined signal in the time window of the present invention are as follows Figure 8 The combined dSNR sequence has a distinct single prominent peak, which can accurately extract the height of the reflecting surface. Compared with the dSNR sequence of a single satellite, the combined dSNR sequence is more conducive to the identification of characteristic frequencies.

[0136] In addition, the different satellite reflection signal fragments used in each time window of the near real-time method come from different azimuths of the water surface, which can more evenly reflect the water surface height; Figure 18 The figure shows the azimuth and elevation angles of the dSNR fragments used in five different time windows of AT01 station DOY201. It can be seen that there are 3 to 6 satellite dSNR fragments in each time window, coming from different azimuths of 0°-220°. Finally, each time window of the near-real-time method contains the ascending and descending tracks of different satellites, which to some extent offsets the effect of the large ascending track inversion value and the small descending track inversion value.

[0137] Comparative analysis of near-real-time spectrum analysis methods and inverse modeling methods:

[0138] The height differences of the reflecting surfaces extracted by the spectrum analysis method and the inverse modeling method in each time window of the C1 and C5 signal groups of the AT01 and SC02 stations are statistically analyzed and plotted on Figure 19 It can be seen that in most time windows, the inversion results of the two methods are very close, and the difference is usually within 0.05m; this shows that both methods have good stability.

[0139] Further statistical analysis revealed that the windows with the largest differences in the inversion results between the two methods are usually those with incomplete time windows in the combined dSNR sequence. In fact, the recombined dSNR sequence within each time window is not always continuous and may have gaps in certain elevation angle ranges. In most cases, these missing data will not affect the correct extraction of the reflector height by the near-real-time spectrum analysis and inverse modeling methods, such as Figure 20As shown; the figure shows the dSNR sequence, LSP analysis and inverse modeling fitting results of the time window combination of AT01 station C5 signal group DOY200 from 5:00 to 5:20; the dSNR sequence with an elevation angle of 9.5°-11.1° in the time window is missing; both the near-real-time spectrum analysis method and the inverse modeling method can be used to invert the water surface height, and the reflector surface heights extracted by the two methods are 13.105m and 13.095m, respectively.

[0140] However, when the dSNR sequence data is missing to a certain extent, the near-real-time spectrum analysis and inverse modeling methods show different inversion performances; Figure 21 The total number of valid windows of the spectrum analysis method and the inverse modeling method during the experiment was counted when the dSNR sequence of the C1 and C5 signal groups of AT01 and SC02 stations was missing more than 1° and more than 2° within the restricted elevation angle range of 5° to 13°. Taking SC02 station as an example, when the dSNR sequence was missing more than 1°, the C1 signal group had a total of 63 valid time windows using the spectrum analysis method, while the inverse modeling method could output 87 valid time windows. For the C5 signal group, under the same conditions, the spectrum analysis method could obtain 58 valid time windows, while the inverse modeling method could obtain 87 valid time windows. The same situation also occurred when the dSNR sequence was missing more than 2°, and the inverse method could obtain more valid time windows than the spectrum analysis method.

[0141] Although LSP can process intermittent and non-uniformly sampled signals, the precision and accuracy of frequency identification will be greatly affected when the signal interruption time is long. In contrast, the inverse modeling method relies on waveform reconstruction. Even if the signal is partially interrupted, as long as the remaining signal waveform is complete, the fitting parameters can still be accurately output. However, when there are too many missing dSNR data in the window, both methods cannot accurately invert the water surface height. In general, the LSP method is more dependent on signal continuity than the inverse modeling method. Therefore, when the available azimuth range of the measuring station is small and there are certain missing points in the combined dSNR sequence within the time window, it is recommended to use the inverse modeling method. When the available azimuth range of the measuring station is large and the combined dSNR sequence within the time window is relatively continuous, the LSP method can obtain more accurate inversion results.

[0142] This paper proposes a GNSS-IR near-real-time water surface height inversion method based on GPS, Galileo, and BDS-3 shared frequency signals. By utilizing the characteristic that dSNR sequences from carriers with the same frequency have similar oscillation waveforms, the method combines the reflection trajectory fragments of each satellite into a complete reflection signal within a small time window, minimizing the inversion time while ensuring high inversion accuracy.

[0143] The sliding step size of the time window determines the temporal resolution of the inversion and can be selected according to the needs of different application scenarios. The study designed a method based on two approaches: spectrum analysis and inverse modeling. Near-real-time water surface height inversion experiments were conducted using observations from the C1 and C5 signal groups at stations AT01 and SC02 to verify the effectiveness of the method.

[0144] The inversion results of near-real-time spectrum analysis and inverse modeling methods are highly consistent with tide observations. AT01 station uses a 20-minute time window and a 5-minute sliding step. The two near-real-time inversion methods have an average window efficiency of 93.4%, an average RMSE of 0.158m, and an average correlation coefficient R of 0.951. SC02 station uses a 30-minute time window and a 10-minute sliding step. The two near-real-time inversion methods have an average window efficiency of 90.1%, an average RMSE of 0.123m, and an average correlation coefficient R of 0.992.

[0145] Compared with the classic GNSS-IR method, the near-real-time spectrum analysis and inverse modeling method can provide more accurate, more real-time, and more uniform temporal resolution water surface height inversion results. Using the near-real-time spectrum analysis and inverse modeling method, the RMSE of the C1, C5, and C1+C5 signal groups at station AT01 decreased by an average of 4.4% compared with the classic method, and the number of inversion points increased by 142.1%. At station SC02, the RMSE of the C1, C5, and C1+C5 signal groups decreased by an average of 18.3%, and the number of inversion points increased by an average of 150.6%.

[0146] Both the near-real-time inverse modeling method and the spectrum analysis method have good stability. When the amount of dSNR sequence data in the time window is small, the near-real-time inverse modeling method can output inversion values in more time windows; when the amount of data in the time window is sufficient, the spectrum analysis method can obtain higher inversion accuracy.

[0147] For stations with an available azimuth range of approximately 180°, a 20-40 minute time window is sufficient for near-real-time solutions, with stable output of high-precision water surface height inversion values over 90% of the time window. With the future modernization of GLONASS, the increase in the number of shared-frequency signals is expected to further shorten the time window for near-real-time inversion methods based on shared-frequency systems, further improving the timeliness of inversions.

[0148] With the above-described preferred embodiments of the present invention as a guide, and with reference to the above description, relevant personnel are fully capable of making various changes and modifications without departing from the technical scope of this invention. The technical scope of this invention is not limited to the contents of the specification and must be determined according to the scope of the claims.

Claims

1. A near real-time water surface height inversion method, characterized in that: The following steps are involved: Step 1: Acquire the shared frequency signal of each navigation system, set different time windows according to different azimuth angles, and obtain the reflected signal sequence fragments of each navigation system at different azimuth angles; Step 2: Use spectrum analysis or inverse modeling to integrate and calculate the fragments of the reflected signal sequence to obtain the height of the reflecting surface, and use the height of the reflecting surface to invert the water surface height.

2. The near real-time water surface height inversion method according to claim 1, characterized in that: Setting different time windows at different azimuths includes: Determine azimuth and elevation ranges; Set the length and step size of the time window; Extract the dSNR of the same carrier frequency of each navigation system; Convert the SNR into linear units and remove the trend term to obtain the dSNR of the reflected signal sequence fragments.

3. The near real-time water surface height inversion method according to claim 1, characterized in that: Inverse modeling methods include: Amplitude normalization is performed on the fragments of the reflected signal sequence; The amplitude-normalized reflected signal fragments are reassembled to obtain: Where, λ is the carrier wavelength; k is the damping coefficient; s is the roughness parameter of the reflecting surface; e is the satellite elevation angle; h is the height of the reflecting surface; C1 and C2 are amplitudes; Perform nonlinear fitting on the amplitude, reflective surface roughness parameters, and reflective surface height to obtain the optimal solution for the reflective surface height; The water surface height is inverted using the optimal reflecting surface height.

4. The near real-time water surface height inversion method according to claim 3, characterized in that: Amplitude normalization uses Hilbert transform to identify the signal envelope and then performs normalization.

5. The near real-time water surface height inversion method according to claim 3, characterized in that: Nonlinear fitting was performed using the nonlinear least squares method.

6. The near real-time water surface height inversion method according to claim 1, characterized in that: Spectral analysis methods include: The dSNR values are averaged within the set elevation angle window to obtain a recombined dSNR sequence; Adaptively denoise the reconstructed dSNR sequence; LSP analysis is used to extract the characteristic frequencies of the recombined dSNR sequence and convert them into water surface height.

7. The near real-time water surface height inversion method according to claim 1, characterized in that: It also includes quality control of the height of the reflective surface.

8. The near real-time water surface height inversion method according to claim 1, characterized in that: Navigation systems include GPS, Galileo and BDS-3.

9. Near real-time water surface height inversion system, characterized by: include: a memory for storing instructions executable by the processor; A processor, configured to execute instructions to implement the near real-time water surface height inversion method according to any one of claims 1 to 8.

10. A computer-readable medium storing computer program code, characterized in that When the computer program code is executed by a processor, the computer program code implements the near real-time water surface height inversion method according to any one of claims 1 to 8.

Images