Real-time inversion method and system for surface water level based on GNSS-R

Abstract

The present application relates to the field of satellite navigation technology, specifically to a GNSS-R based real-time ground water level inversion method and system, comprising the following steps: collecting multi-band GNSS signals, combining an elevation model to solve the specular reflection point, extracting the delay Doppler characteristics of the reflected signal, using environmental parameters to compensate for the calculation of the vertical distance from the antenna to the water surface, and combining the antenna geodetic height to generate the absolute water level of the ground water body. In the present application, by synchronously collecting multi-constellation data and achieving time alignment, the data fusion accuracy and response speed are improved, the clock error is corrected by analyzing the orbit parameters, the reflection points are accurately selected by combining the elevation model to exclude interference, the delay Doppler grid is constructed to perform coherent integration, the peak value characteristics are extracted to enhance the spectrum resolution and capture robustness, the height error is compensated by introducing multi-source environmental parameters, the geodetic datum is applied to unify the reference, and the accuracy and continuous stability of the centimeter-level absolute water level monitoring are ensured, effectively supporting real-time application of large-scale water bodies.

Description

Technical Field

[0001] This invention relates to the field of satellite navigation technology, and in particular to a method and system for real-time inversion of surface water level based on GNSS-R. Background Technology

[0002] With the intensification of global climate change and human activities, the rational management and protection of water resources have become particularly important. Accurate and real-time monitoring of surface water changes is crucial for preventing natural disasters such as floods and droughts, and for the scientific planning of water resource utilization. However, traditional methods such as buoy methods and pressure sensor methods, while providing precise data, face challenges such as high installation and maintenance costs (annual maintenance costs may exceed tens of thousands of US dollars), limited coverage (usually limited to specific locations), and inability to achieve large-scale continuous monitoring. Surface water level inversion equipment based on Global Navigation Satellite System Reflectance (GNSS-R) technology has emerged to address this need, obtaining water surface height information by receiving and analyzing the reflected waves of GNSS signals on the Earth's surface.

[0003] This technology boasts all-weather, all-time operation capabilities, high spatiotemporal resolution, and low cost advantages, making it particularly suitable for continuous observation of large-scale water bodies in complex geographical environments. With the continuous improvement of GNSS systems and the rapid development of related fields such as microwave electronics and digital signal processing, the widespread availability of low-cost GNSS receivers and high-performance computing platforms has made it possible to construct a complete surface water level inversion system. This equipment aims to meet diverse needs from scientific research to government management and even private enterprises, demonstrating broad application prospects in multiple fields such as scientific research, disaster early warning, resource management, and environmental protection. It also adopts a modular design to adapt to different user requirements.

[0004] Therefore, there is an urgent need to develop a low-cost, real-time, and high-precision surface water level inversion method and system based on GNSS-R to overcome the shortcomings of traditional methods. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of existing technologies and propose a real-time surface water level inversion method based on GNSS-R.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: a real-time surface water level inversion method based on GNSS-R, comprising the following steps:

[0007] S1: Detect direct and reflected signals from L-band GNSS satellites, simultaneously collect signal-to-noise ratio observations, carrier phase observations, and satellite pseudorange data under multi-constellation joint observations, and generate a multi-band GNSS raw observation signal set;

[0008] S2: Based on the analysis of the transmission time of the original observation signal set of the multi-band GNSS, the orbital parameters are called to calculate the geometric transmission path, and the theoretical reflection position is screened using the digital elevation model and water body mask to obtain the coordinates of the water surface mirror reflection point;

[0009] S3: For the region pointed to by the coordinates of the water surface mirror reflection point, calculate the correlation value between the reflection component and the copy code of the multi-band GNSS raw observation signal set, perform coherent integration and extract the peak delay, and obtain the delay Doppler map features of the reflection signal;

[0010] S4: Call the delayed Doppler image feature of the reflected signal to read the time difference value, combine the coordinates of the water surface mirror reflection point to calculate the path difference and elevation angle sine value, use atmospheric, tidal and wind wave parameters to compensate for the foundation height, and establish the vertical distance value from the antenna phase center to the water surface.

[0011] S5: Obtain the geodetic elevation parameters of the antenna reference ellipsoid based on the vertical distance from the antenna phase center to the water surface, calculate the difference between the geodetic elevation parameters of the reference ellipsoid and the vertical distance value, and generate the absolute water level elevation of the surface water body.

[0012] As a further aspect of the present invention, the multi-band GNSS raw observation signal set includes a signal-to-noise ratio observation sequence, a carrier phase observation sequence, and multi-constellation pseudorange observation data;

[0013] The coordinates of the water surface mirror reflection point include the geodetic coordinates of the mirror reflection point, the validity marker of the reflection point, and the boundary data of the water body mask.

[0014] The delayed Doppler graph features of the reflected signal include the time delay at the specular reflection point, the Doppler frequency shift at the specular reflection point, and the peak power of the delayed Doppler graph;

[0015] The vertical distance from the antenna phase center to the water surface includes the two-way propagation path difference, atmospheric and geophysical corrections, and antenna height inversion value.

[0016] The absolute water level elevation of the surface water body includes the station ellipsoid height, the absolute water surface elevation, and the water level time series.

[0017] As a further aspect of the present invention, the satellite pseudorange data is a distance metric obtained by multiplying the signal propagation time by the speed of light, and is used for timing alignment of the signal reception time.

[0018] Digital elevation models are data containing surface elevation information, used to perform elevation anomaly correction on theoretical specular reflection points;

[0019] The water mask is used to distinguish between land and water geographical data, determine the reflection point landing area attributes, and remove invalid areas.

[0020] As a further aspect of the present invention, the peak delay is the maximum power point corresponding quantity extracted from the two-dimensional delayed Doppler spectrum, which is used to characterize the time interval between the arrival of the reflected and direct signals.

[0021] The vertical distance value is the vertical height from the antenna phase center to the water surface, obtained after time interval conversion, geometric path correction, and environmental parameter compensation.

[0022] The reference ellipsoid geodetic elevation parameters are obtained by calling the receiver's precise single-point positioning results and filtering and smoothing them to determine the absolute height of the station antenna phase center relative to the reference ellipsoid surface.

[0023] As a further aspect of the present invention, the step of obtaining S1 is as follows:

[0024] The system detects right-hand circularly polarized direct signals emitted by L-band GNSS satellites and left-hand circularly polarized reflected signals reflected by the water surface, separates electromagnetic wave energy in different polarization directions, locks the signal frequency, and obtains the basic signal flow of the dual-polarization channel.

[0025] The carrier and ranging code information of the basic signal stream of the dual-polarized channel are called to demodulate the signal-to-noise ratio (SNR) and carrier phase observations of the direct and reflected links. The validity of the SNR and carrier phase observations is screened to obtain the observation characteristic parameters of the direct and reflected links.

[0026] Based on the observation characteristic parameters of the direct and reflected links, satellite pseudorange data under multi-constellation joint observation are collected synchronously. The satellite pseudorange data and characteristic parameters are time-aligned according to the reception time to establish a multi-band GNSS raw observation signal set.

[0027] As a further aspect of the present invention, the step of obtaining S2 is as follows:

[0028] Based on the multi-band GNSS raw observation signal set, the satellite signal transmission time is analyzed, the precise orbital parameters and satellite clock error parameters are called to correct the satellite clock deviation, the three-dimensional orbital position of the satellite's center of mass at the transmission time is calculated, and the instantaneous spatial position parameters of the satellite are obtained.

[0029] The instantaneous spatial position parameters of the satellite and the preset geodetic coordinate parameters of the receiver are used to construct a signal transmission geometric link. The theoretical position where the signal propagation path is tangent to the reference ellipsoid is calculated based on the principle of geometric optics reflection, and the spatial vector of the theoretical specular reflection point is generated.

[0030] The spatial vector of the theoretical mirror reflection point is used to introduce digital elevation model data for elevation anomaly correction. Water body mask data is used to determine the surface attributes of the reflection point's landing area and remove invalid areas to establish the coordinates of the water surface mirror reflection point.

[0031] As a further aspect of the present invention, the step of obtaining S3 is as follows:

[0032] For the area pointed to by the coordinates of the water surface mirror reflection point, the reflected signal data segment is extracted from the multi-band GNSS raw observation signal set, a local signal copy with the same pseudo-random code structure as the satellite transmission signal is generated, and the reflected signal components and the local copy code sequence are obtained.

[0033] The reflected signal component and the local copy code sequence are called to construct a two-dimensional search grid in the time delay domain and Doppler frequency domain. Cross-correlation calculation is performed on the reflected signal component and the local copy code sequence and coherence integral is accumulated to obtain the signal power response intensity of each grid node and generate a two-dimensional delayed Doppler power spectrum array.

[0034] Based on the two-dimensional delayed Doppler power spectrum array, the power values ​​of each grid node are traversed, the index of the location of the point with the maximum power is retrieved, and the time delay parameter is extracted as a feature quantity to obtain the delayed Doppler map feature of the reflected signal.

[0035] As a further aspect of the present invention, the step of obtaining S4 is as follows:

[0036] The delayed Doppler map features of the reflected signal are called up, the time difference between the direct signal and the reflected signal arriving at the receiver is read, and a geometric link model is constructed by combining the coordinates of the water surface mirror reflection point and the satellite spatial coordinates. The geometric difference of the two-way path and the sine of the incident elevation angle of the satellite relative to the reflection plane are calculated to obtain the dual-station geometric delay and angle dataset.

[0037] The time difference is converted into a distance by calling the bi-station geometric delay and angle dataset and introducing the vacuum speed of light constant. The geometric height is calculated by using the bi-path geometric difference and the sine of the incident elevation angle. The geometric path is corrected by combining the Earth curvature parameter to generate the geometric inversion foundation vertical height.

[0038] Based on the geometric inversion foundation vertical height, atmospheric delay environment calibration parameters, tidal environment influence parameters, and wind and wave environment influence parameters are obtained. The error superposition and cancellation value of the foundation vertical height is calculated and dynamically compensated, and the vertical distance value from the antenna phase center to the water surface is established.

[0039] As a further aspect of the present invention, the step of obtaining S5 is as follows:

[0040] The observation epoch is determined based on the vertical distance from the antenna phase center to the water surface. The sequence of precise single-point positioning results from the receiver is retrieved. Points with jumps in positioning coordinates are removed by filtering and smoothing, and the average vertical height is calculated to obtain the reference geodetic elevation parameters of the station.

[0041] The difference between the station's reference geodetic elevation parameters and the vertical distance from the antenna phase center to the water surface is calculated by calling the station's reference geodetic elevation parameters and the vertical distance from the antenna phase center to the water surface. The attitude error of the difference is corrected by the equipment installation tilt angle and signal quality factor, and a weighted calculation is performed to generate a dynamic elevation settlement value for the water surface.

[0042] Based on the calculated dynamic elevation values ​​of the water surface, the parameters of the geoid refinement model for this region are introduced to calculate the elevation anomaly separation between the reference ellipsoid and the geoid. The water surface elevation is then converted to an orthographic system and unified with the benchmark to establish the absolute water level elevation of the surface water body.

[0043] The GNSS-R-based real-time surface water level inversion system includes: a sensing end signal acquisition module, a mirror point spatial positioning module, a signal feature calculation module, a vertical ranging inversion module, and an absolute water level calculation module.

[0044] The sensing end signal acquisition module is used to capture L-band signals through a dual-polarized antenna, acquire signal-to-noise ratio and carrier phase data through a radio frequency front-end and a digital receiver, synchronize multiple constellation pseudoranges, and establish a multi-band GNSS raw observation signal set.

[0045] The mirror point spatial positioning module is used to call the original observation signal set of the multi-band GNSS to analyze the transmission time, calculate the geometric path in combination with the precise orbit parameters, use the digital elevation model to screen the landing area, and generate the coordinates of the water surface mirror reflection point.

[0046] The signal feature calculation module is used to drive the edge computing unit to call the multi-band GNSS raw observation signal set to extract data for the area pointed to by the coordinates of the water surface mirror reflection point, generate a local copy code, perform cross-correlation calculation in the time delay and Doppler domain and accumulate coherent integrals to construct a two-dimensional power spectrum, traverse the grid nodes to retrieve the point with the maximum power, extract the time delay parameters, and obtain the delayed Doppler map features of the reflected signal.

[0047] The vertical ranging inversion module is used to read the time difference by utilizing the delayed Doppler image features of the reflected signal, calculate the geometric path difference by combining the coordinates of the water surface mirror reflection point, introduce environmental parameters to compensate for the foundation height, and generate the vertical distance value from the antenna phase center to the water surface.

[0048] The absolute water level calculation module is used to retrieve locally stored station reference parameters based on the vertical distance from the antenna phase center to the water surface, calculate the elevation difference, introduce a unified benchmark from the geoid model, and establish the absolute water level elevation of the surface water body.

[0049] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0050] In this invention, multi-constellation signal-to-noise ratio observations, carrier phase observations, and pseudorange data are simultaneously acquired to form a multi-band signal set, achieving polarization separation and time alignment, improving data fusion accuracy and response speed. Precise orbital parameters are used to correct clock errors by analyzing the launch time, calculating satellite positions to construct geometric links, and using a digital elevation model combined with a water mask to screen reflection points, enhancing position accuracy and eliminating interference. A local copy code sequence is generated for the reflection area, and a time-delay Doppler grid is constructed to perform cross-correlation, accumulation, and coherent integration, extracting peak time delay frequency shift features to improve spectral resolution and acquisition robustness. The time difference is read to calculate the sine value of the two-way path difference elevation angle, and atmospheric tidal and wave parameters are introduced for compensation to reduce the impact of environmental errors on altitude estimation. Based on the calculated vertical distance difference of the reference ellipsoid elevation, an attitude-corrected geoid model is applied to unify the benchmark, ensuring centimeter-level absolute water level accuracy and continuous monitoring stability, supporting applications in large-scale rivers, lakes, and urban water bodies. Attached Figure Description

[0051] Figure 1 This is a flowchart of the main steps of the present invention;

[0052] Figure 2 This is a schematic diagram of the mirror reflection point positioning of the present invention;

[0053] Figure 3 This is a time series comparison chart of the inverted water level and the reference value according to the present invention;

[0054] Figure 4 This is a scatter plot of the inversion results and measured data of the present invention;

[0055] Figure 5 This is a comparison chart of the high-resolution tidal inversion of the present invention and NOAA's actual measurements;

[0056] Figure 6 This is a scatter plot of GPS inversion and NOAA tide levels from the present invention.

[0057] Figure 7 This is a schematic diagram of the two-dimensional delayed Doppler power spectrum of the present invention;

[0058] Figure 8 This is a schematic diagram of a reservoir scene according to the present invention;

[0059] Figure 9 This is a schematic diagram of a river scene according to the present invention;

[0060] Figure 10 This is a schematic diagram of an urban flooding scenario according to the present invention;

[0061] Figure 11 This is a hardware structure diagram of the present invention;

[0062] Figure 12 This is a system flowchart of the present invention. Detailed Implementation

[0063] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0064] Please see Figures 1-10 A real-time surface water level inversion method based on GNSS-R includes the following steps:

[0065] S1: Detect right-hand circularly polarized direct signals emitted by L-band GNSS satellites and left-hand circularly polarized reflected signals reflected by the water surface. Simultaneously acquire signal-to-noise ratio observations, carrier phase observations, and satellite pseudorange data under multi-constellation joint observations of the direct and reflected signals to generate a multi-band GNSS raw observation signal set.

[0066] S2: Based on the analysis of satellite signal transmission time from the multi-band GNSS raw observation signal set, the precise orbit parameters and satellite clock error parameters are called, and the geometric transmission path of the satellite signal is calculated by combining the geodetic coordinate parameters preset by the receiver. The theoretical reflection position is spatially filtered using digital elevation model data and water mask data to generate the coordinates of the water surface mirror reflection point.

[0067] S3: For the region pointed to by the coordinates of the water surface mirror reflection point, calculate the correlation value between the reflected signal component and the local copy code in the multi-band GNSS raw observation signal set, perform two-dimensional coherent integration in the time delay domain and Doppler frequency domain, extract the time delay feature corresponding to the point with the maximum power in the spectrum, and obtain the delayed Doppler map feature of the reflected signal.

[0068] S4: Call the delayed Doppler image feature of the reflected signal to read the time difference between the direct signal and the reflected signal reaching the phase center of the receiver antenna. Combine the coordinates of the water surface mirror reflection point and the satellite spatial coordinates to construct a geometric relationship model. Calculate the geometric difference of the two-way path and the sine value of the incident elevation angle of the satellite relative to the reflection plane. Obtain atmospheric delay environment calibration parameters, tidal environment influence parameters, and wind and wave environment influence parameters. Calculate the vertical height of the foundation based on the geometric difference of the two-way path and the sine value of the incident elevation angle. Use the atmospheric delay environment calibration parameters, tidal environment influence parameters, and wind and wave environment influence parameters to dynamically compensate the vertical height of the foundation and generate the vertical distance value from the antenna phase center to the water surface.

[0069] S5: Based on the vertical distance from the antenna phase center to the water surface, obtain the geodetic elevation parameters of the receiver antenna phase center reference ellipsoid, calculate the difference between the geodetic elevation parameters of the receiver antenna phase center reference ellipsoid and the vertical distance from the antenna phase center to the water surface, and establish the absolute water level elevation of the surface water body.

[0070] The multi-band GNSS raw observation signal set includes signal-to-noise ratio observation sequence, carrier phase observation sequence, and multi-constellation pseudorange observation data;

[0071] The coordinates of the water surface mirror reflection point include the geodetic coordinates of the mirror reflection point, the validity marker of the reflection point, and the boundary data of the water body mask.

[0072] The characteristics of the delayed Doppler graph of the reflected signal include the time delay at the specular reflection point, the Doppler frequency shift at the specular reflection point, and the peak power of the delayed Doppler graph;

[0073] The vertical distance from the antenna phase center to the water surface includes the two-way propagation path difference, atmospheric and geophysical corrections, and antenna height inversion value.

[0074] The absolute water level elevation of a surface water body includes the station ellipsoid height, the absolute water surface elevation, and the water level time series;

[0075] Please see Figure 1 The steps to obtain S1 are as follows:

[0076] The system detects right-hand circularly polarized direct signals emitted by L-band GNSS satellites and left-hand circularly polarized reflected signals reflected by the water surface, separates electromagnetic wave energy in different polarization directions, locks the signal frequency, and obtains the basic signal flow of the dual-polarization channel.

[0077] First, using a dual-polarized antenna array at the receiving end, the system locks onto the right-hand circularly polarized direct wave emitted by an L-band GNSS satellite and the left-hand circularly polarized wave reflected from the water surface, targeting the zenith and ground reflection directions respectively. The front-end center frequency is configured as 1575.42MHz to correspond to the L1 band or 1176.45MHz to correspond to the L5 band. The received analog signal is processed by a low-noise amplifier with the gain controlled between 30dB and 50dB. Then, it is down-converted to an intermediate frequency using a quadrature demodulation circuit, and analog-to-digital conversion is performed at a rate more than four times the intermediate frequency, such as 20MHz, to quantize the in-phase and quadrature components. The power level of each channel is monitored in real time, and a background noise level baseline is set. The standard value is the average thermal noise of the receiver in the absence of a signal, for example, -130dBm. Based on the signal detection probability curve, the direct signal determination coefficient is set to 10dB, and the direct signal acquisition threshold is calculated as -130dBm plus 10dB, which is -120dBm. Similarly, the reflection signal determination coefficient is set to 3dB, and the reflection acquisition threshold is calculated as -127dBm. When the level of the direct channel is higher than -120dBm and the level of the reflection channel is higher than -127dBm, the local oscillator is triggered to lock the carrier frequency and keep the deviation within ±500Hz. The digital streams carrying carrier and ranging code information are output respectively, thereby obtaining the basic signal stream of the dual-polarized channel.

[0078] By calling the basic signal stream of the dual-polarized channel to demodulate the carrier and ranging code information, the signal-to-noise ratio (SNR) and carrier phase observations of the direct and reflected links are extracted. The numerical validity of the SNR and carrier phase observations is then screened to obtain the observation characteristic parameters of the direct and reflected links.

[0079] First, cross-correlation is performed on the basic signal stream of the dual-polarized channel using a locally generated pseudo-random noise code sequence. At integration durations of 1 ms or 20 ms, the sum of squares of the integral amplitudes of the in-phase and quadrature branches is calculated. Then, the carrier power-to-noise ratio (SNR), i.e., the signal-to-noise ratio, is calculated based on the noise floor. Simultaneously, the phase register of the numerically controlled oscillator within the carrier tracking loop is read. Carrier phase observations are constructed by splicing full-cycle counts and fractional-period phases. During the data filtering stage, based on receiver sensitivity statistics, the direct signal-to-noise ratio reference value is set to 35 dB-Hz with a 3 dB safety margin. The calculated... The effective threshold for the signal-to-noise ratio (SNR) of the direct signal is 38 dB-Hz. Similarly, based on the attenuation characteristics of the reflected signal, the effective threshold for the SNR of the reflected signal is set to 25 dB-Hz. The measured direct link value of 45 dB-Hz is compared with the threshold and determined to be valid. The reflected link value of 28 dB-Hz, which is higher than the threshold, is retained. For the carrier phase, the cycle slip judgment benchmark value is calculated based on the wavelength of the L1 band, which is about 19 cm. The phase difference threshold is set to half the wavelength, i.e., 9.5 cm. If the phase difference conversion distance between consecutive epochs exceeds 9.5 cm, it is regarded as a cycle slip and discarded. The observation characteristic parameters of the direct and reflected links are obtained.

[0080] Based on the observation characteristic parameters of direct and reflected links, satellite pseudorange data under multi-constellation joint observation are collected synchronously. The satellite pseudorange data and characteristic parameters are time-aligned according to the reception time to establish a multi-band GNSS raw observation signal set.

[0081] First, based on the observation characteristic parameters of direct and reflected light links, the multi-channel baseband processing unit is driven to analyze the navigation messages of the BeiDou, GPS, and Galileo systems. The signal propagation time is calculated by combining the receiver's local time and the transmission time of each satellite, and multiplied by the speed of light in vacuum, 299,792,458 m / s, to generate satellite pseudorange data. The time synchronization tolerance window is set to 1 microsecond based on the stability parameters of the receiver's clock crystal oscillator. The receiving time tag in the characteristic parameters is read, for example, the second within the week, 123,456.789,000. Satellite pseudorange data streams with time differences within this range are traversed and retrieved. For example, if the time tag of a pseudorange data differs from the characteristic parameter by 0.5 microseconds, which is less than the 1 microsecond threshold, a successful match is determined and the data lines are merged. If it exceeds the window, interpolation is performed to achieve time alignment. After removing isolated points, the signal-to-noise ratio, carrier phase, and pseudorange data are reconstructed according to the satellite number and frequency point to establish a multi-band GNSS raw observation signal set.

[0082] Please see Figures 1-2The steps to obtain S2 are as follows:

[0083] Based on the analysis of satellite signal transmission time using multi-band GNSS raw observation signal set, the satellite clock deviation is corrected by calling precise orbital parameters and satellite clock error parameters, the three-dimensional orbital position of the satellite's center of mass at the transmission time is calculated, and the instantaneous spatial position parameters of the satellite are obtained.

[0084] First, based on the receiver time stamp in the observation data packet, such as the second of this week 345600.0, the corresponding satellite clock error coefficient, such as 120 microseconds, is extracted from the precise clock error file published by IGS. The precise signal transmission time is analyzed by subtracting the pseudorange propagation time and clock error correction value from the reception time. Then, the discrete coordinates in the precise ephemeris file are called, and polynomial interpolation is performed on 5 epochs before and after the transmission time to calculate the three-dimensional coordinates in the Earth-fixed coordinate system, such as X-axis 25000000.5 meters, Y-axis 10000000.2 meters, and Z-axis 5000000.1 meters. The Earth's rotation angular velocity of 7.292e-5 radians / second is introduced and combined with the propagation time to construct a rotation matrix to complete the coordinate rotation transformation to eliminate the rotation effect and obtain the instantaneous spatial position parameters of the satellite.

[0085] The signal transmission geometric link is constructed by calling the satellite's instantaneous spatial position parameters and the receiver's preset geodetic coordinate parameters. The theoretical position where the signal propagation path is tangent to the reference ellipsoid is calculated based on the principle of geometric optics reflection, and the spatial vector of the theoretical mirror reflection point is generated.

[0086] First, using the semi-major axis of the WGS-84 reference ellipsoid at 6378137 meters as a reference, a geometric link endpoint is constructed using the satellite's instantaneous spatial position parameters and the preset receiver geodetic coordinates. Based on the geometric optics principle that the incident angle equals the reflection angle, a set of nonlinear equations for the incident and reflection vectors is established. A differential iterative algorithm is used to search for points on the reference ellipsoid that satisfy the conditions of coplanarity and equiangularity. Based on the positioning accuracy requirements, the iteration convergence threshold is set to 0.001 meters. During the iterative calculation, the Euclidean distance between the coordinates of the current iteration point and the coordinates of the previous iteration point is calculated. For example, if the current difference is 0.0005 meters, which is less than the 0.001-meter threshold, the result is determined to be converged and the iteration is terminated. The XYZ components of the point in the Earth-fixed coordinate system are output, generating the spatial vector of the theoretical specular reflection point.

[0087] To correct elevation anomalies, the spatial vector of the theoretical mirror reflection point is introduced into digital elevation model data. Water body mask data is used to determine the surface properties of the reflection point's landing area and remove invalid areas to establish the coordinates of the water surface mirror reflection point.

[0088] First, the spatial vector of the theoretical mirror reflection point is converted into geodetic latitude and longitude, such as 32.5 degrees north latitude and 118.8 degrees east longitude. Then, a digital elevation model with a spatial resolution of 30 meters is introduced for elevation anomaly correction. The surface elevation at this coordinate is read, for example, 15.6 meters. The reflection point is translated to this elevation surface along the normal direction. The geometric path intersections from the satellite to the surface and from the surface to the receiver are recalculated. Then, the binarized water mask database is called to retrieve the raster pixel values ​​where the corrected coordinates fall. The effective water area interpretation benchmark value is set to 255. If the retrieved pixel value is 0, i.e., land, it is discarded. If it is 255, which is equal to the benchmark value, it is determined to be an effective water area and retained. In this way, the coordinates of the water surface mirror reflection point are established.

[0089] Please see Figure 1 , Figure 7 The steps to obtain S3 are as follows:

[0090] For the region pointed to by the coordinates of the water surface mirror reflection point, the reflected signal data segment is extracted from the multi-band GNSS raw observation signal set, a local signal copy with the same pseudo-random code structure as the satellite transmitted signal is generated, and the reflected signal components and the local copy code sequence are obtained.

[0091] First, for the area pointed to by the coordinates of the water surface reflection point, calculate the length difference between the geometric path of the reflected signal and the direct path, for example, 150 meters. Divide this by the speed of light in vacuum to obtain a theoretical delay of about 500 nanoseconds. Based on the signal symbol period, set the search window width coefficient to 1.0. Calculate the search window width as 1023 chips multiplied by the coefficient, i.e., 1023 code width. Extract a 20-millisecond intermediate frequency digital sampling segment of the reflection channel from the multi-band GNSS raw observation signal set. Identify the current visible satellite pseudo-random code, such as PRN12. Configure the local signal generator to generate a standard pseudo-random noise code sequence of 1023 chips with a clock frequency of 1.023MHz. Adjust its starting phase to match the search window. Output the discrete sampling data of the reflected signal and the local copy to obtain the reflected signal components and the local copy code sequence.

[0092] By calling the reflected signal component and the local copy code sequence, a two-dimensional search grid in the time delay domain and Doppler frequency domain is constructed. Cross-correlation calculation is performed on the reflected signal component and the local copy code sequence, and coherence integral is accumulated to obtain the signal power response intensity of each grid node, thereby generating a two-dimensional delayed Doppler power spectrum array.

[0093] First, with the theoretical delay set to zero, a two-dimensional search grid is constructed with 21 nodes in the delay range of -1.0 to +1.0 chip increments of 0.1 chip increments and 21 nodes in the Doppler frequency shift range of -500Hz to +500Hz increments of 50Hz. The reflected signal component and the local copy code sequence are called. For each node, the local code is cyclically shifted and the carrier phase rotation factor corresponding to the frequency is superimposed. Pointwise complex multiplication is performed with the reflected signal component. The product results are accumulated and summed with a 1-millisecond integration time. The power response intensity is calculated by the square of the modulus, for example, 4500 amplitude units. This is then filled into the corresponding cell to generate a two-dimensional delayed Doppler power spectrum array.

[0094] Based on the two-dimensional delayed Doppler power spectrum array, the power values ​​of each grid node are traversed, the index of the location of the point with the maximum power is retrieved, and the time delay parameter is extracted as a feature quantity to obtain the delayed Doppler map feature of the reflected signal.

[0095] First, initialize the peak power register to 0. Traverse the power data of 441 grid nodes in the two-dimensional delayed Doppler power spectrum array. Compare the current node value, for example, 4800, with the register value. If 4800 is greater than 0 in the register, update the register and record the time delay axis index, such as 11, and the frequency axis index, such as 10, until the traversal ends and the position of the maximum power point is locked. Map this position to a time delay parameter, such as 0.0 chip, as a key quantity characterizing the arrival time difference of the reflected signal relative to the direct signal, and obtain the delayed Doppler map features of the reflected signal.

[0096] Please see Figure 1 , Figure 7 The steps to obtain S4 are as follows:

[0097] The delayed Doppler image features of the reflected signal are called up, the time difference between the direct signal and the reflected signal arriving at the receiver is read, and a geometric link model is constructed by combining the coordinates of the water surface mirror reflection point and the satellite spatial coordinates. The geometric difference of the two-way path and the sine of the incident elevation angle of the satellite relative to the reflection plane are calculated to obtain the dual-station geometric delay and angle dataset.

[0098] First, read the time delay value corresponding to the relevant peak in the delayed Doppler image feature of the reflected signal, such as 10.5 chips, and convert it into a time difference value. Combine the instantaneous spatial coordinates of the satellite and the coordinates of the water surface reflection point, calculate the Euclidean distances from the satellite to the reflection point (24,000 km), from the reflection point to the receiver (25,000 km), and from the satellite to the receiver (20 m). Obtain the geometric difference of the two-way path through addition and subtraction operations. At the same time, calculate the cosine of the angle between the satellite position vector and the surface normal vector as the sine of the incident elevation angle, such as 0.5. Package the above parameters to obtain the bistatic geometric delay and angle dataset.

[0099] By calling the bi-station geometric delay and angle dataset and introducing the vacuum speed of light constant, the time difference is converted into a distance quantity. The geometric height is calculated using the bi-path geometric difference and the sine of the incident elevation angle. The geometric path is corrected by combining the Earth curvature parameter, and the geometric inversion foundation vertical height is generated.

[0100] First, the time difference is converted into a distance quantity by introducing the vacuum speed of light constant from the bistatic geometric delay and angle dataset, using the formula:

[0101] ,

[0102] in, The vertical height of the geometric inversion foundation is in meters. This is the time delay, expressed in seconds. Represents the speed of light in vacuum, measured in meters per second. This is the equivalent distance correction for system hardware delay, in meters. This is the quadratic term of the geometric correction for the curvature of the Earth, in square meters; The value of the sine of the incident elevation angle is a dimensionless parameter. The time constant for low elevation angle multipath deviation is expressed in seconds. This represents the number of sampling points for the fitted residuals. For sampling point index; The residuals for fitting the geometric model are in meters. This is the residual weighting factor, i.e., a dimensionless parameter;

[0103] Formula calculation process: First, the equivalent distance correction amount of system hardware delay is used. After eliminating inherent equipment biases in the observation distance, and incorporating the quadratic term of the geometric correction for Earth's curvature... Geometric correction is performed on long-distance transmission paths, followed by utilization of the low elevation angle multipath bias time constant. With observation time delay The ratio is used to construct a dimensionless correction factor, which is used to fine-tune the incident geometry term in the denominator. Finally, by subtracting the weighted geometric model fitting residual term, the preliminary inversion height containing only geometric path information is obtained, and the geometric inversion base vertical height is generated.

[0104] Based on the geometric inversion of the basic vertical height, atmospheric delay environment calibration parameters, tidal environment influence parameters, and wind and wave environment influence parameters are obtained. The error superposition and cancellation value of the basic vertical height is calculated and dynamic compensation is performed. The vertical distance value from the antenna phase center to the water surface is established.

[0105] First, the station location is determined based on the geometric inversion of the basic vertical height. The Saastamoinen model is called, and meteorological parameters, namely air pressure 1013 hPa and air temperature 20℃, are input. The tropospheric delay calibration amount is calculated, for example, 0.05 meters. The vertical displacement parameter of 0.02 meters is obtained by querying the FES2014 global tidal model. At the same time, the wind and wave deviation calculation coefficient is set to 0.001 according to the sea surface roughness model. Real-time wind speed data of 10 m / s is read, and the wind and wave roughness deviation factor is calculated as wind speed squared × coefficient = 100 × 0.001 = 0.1 meters. The basic vertical height (assuming the above formula calculation result is 15.5 meters) is subjected to addition and subtraction operations, that is, adding delay and subtracting wind and wave deviation from tide, to obtain a compensated result of 15.47 meters. The vertical distance value from the antenna phase center to the water surface is established.

[0106] Please see Figure 1 , Figures 3-10 The steps to obtain S5 are as follows:

[0107] The observation epoch is determined based on the vertical distance from the antenna phase center to the water surface. The sequence of precise single-point positioning results from the receiver is retrieved. Points with jumps in positioning coordinates are removed by filtering and smoothing, and the average vertical height is calculated to obtain the reference geodetic elevation parameters of the station.

[0108] First, using the observation epoch associated with the vertical distance from the antenna phase center to the water surface, such as 345600 seconds within a week, as the anchor point, the precise single-point positioning sequence of the previous hour in the receiver memory is traced back. Based on the statistical characteristics of the station's historical positioning data, an elevation jump detection threshold is set. For example, if the historical elevation standard deviation is 0.16 meters, 0.16 × 3 = 0.48 meters is calculated according to the principle of 3 times the standard deviation. The threshold is then rounded down to 0.5 meters. The elevation difference between adjacent epochs is calculated by traversing the sequence. If the difference of a certain epoch is 0.8 meters, which is greater than the 0.5-meter threshold, it is discarded. The arithmetic mean is then performed on the vertical coordinate components of the remaining valid data to calculate the average geodetic elevation under the WGS-84 system, for example, 50.0 meters, thus obtaining the station's reference geodetic elevation parameters.

[0109] The station reference geodetic elevation parameters and the vertical distance from the antenna phase center to the water surface are called up. The difference between the station reference geodetic elevation parameters and the vertical distance from the antenna phase center to the water surface is calculated. The attitude error of the difference is corrected by the equipment installation tilt angle and signal quality factor and weighted calculation is performed to generate the dynamic elevation settlement value of the water surface.

[0110] First, a subtraction operation is performed between the station's reference geodetic elevation parameter of 50.0 meters and the vertical distance from the antenna phase center to the water surface of 15.56 meters, yielding a preliminary water surface ellipsoid elevation of 34.44 meters. Then, the antenna installation tilt angle measured by an external attitude sensor is introduced, for example, 2.5 degrees. The vertical distance is geometrically corrected to 15.576 meters using the cosine value of 0.999. The weighted attenuation calculation formula is set as follows:

[0111] ,

[0112] Among them, the ideal specular reflection reference value Set to 45dB, tolerance range Set to 75dB, substitute the observed signal-to-noise ratio of 30dB to calculate the weighting coefficient W=1−∣45−30∣ / 75=0.8, apply this weight to the corrected difference for weighted smoothing, and generate the dynamic elevation settlement value of the water surface.

[0113] Based on the numerical calculation of dynamic water surface elevation, the parameters of the geoid refinement model in this region are introduced to calculate the elevation anomaly separation between the reference ellipsoid and the geoid. The water surface elevation is then converted to an orthographic system and unified with the benchmark to establish the absolute water level elevation of the surface water body.

[0114] First, based on the calculated dynamic elevation of the water surface at 34.38 meters, a refined geoid model for the local area is loaded. According to the 114.3-degree latitude and longitude and 30.5-degree index grid of the station, the elevation anomaly separation between the reference ellipsoid and the geoid is extracted, for example, -0.42 meters. Following the principle of orthographic height calculation, the anomaly is subtracted from the calculated value, i.e., 34.38 meters minus -0.42 meters, to obtain 34.80 meters. This completes the benchmark conversion from geometric ellipsoidal height to physical orthographic height system, and establishes the absolute water level elevation of the surface water body.

[0115] Please see Figures 11-12 A GNSS-R-based real-time surface water level inversion system includes:

[0116] The sensor terminal signal acquisition module, the mirror point spatial positioning module, the signal feature calculation module, the vertical ranging inversion module, and the absolute water level calculation module are all included.

[0117] The sensing end signal acquisition module is used to capture L-band signals through a dual-polarized antenna, acquire signal-to-noise ratio and carrier phase data through a radio frequency front-end and a digital receiver, synchronize multiple constellation pseudoranges, and establish a multi-band GNSS raw observation signal set.

[0118] The mirror point spatial positioning module is used to call the multi-band GNSS raw observation signal set to analyze the transmission time, combine it with precise orbit parameters to calculate the geometric path, use the digital elevation model to screen the landing area, and generate the coordinates of the water surface mirror reflection point.

[0119] The signal feature calculation module is used to drive the edge computing unit to call the multi-band GNSS raw observation signal set to extract data for the region pointed to by the coordinates of the water surface mirror reflection point, generate local copy code, perform cross-correlation calculation in the time delay and Doppler domain and accumulate coherent integrals to construct a two-dimensional power spectrum, traverse the grid nodes to search for the point with the maximum power, extract the time delay parameters, and obtain the delayed Doppler map features of the reflected signal.

[0120] The vertical ranging inversion module is used to read the time difference of the delayed Doppler image features of the reflected signal, calculate the geometric path difference by combining the coordinates of the water surface mirror reflection point, introduce environmental parameters to compensate for the foundation height, and generate the vertical distance value from the antenna phase center to the water surface.

[0121] The absolute water level calculation module is used to retrieve locally stored station reference parameters based on the vertical distance from the antenna phase center to the water surface, calculate the elevation difference, introduce the geoid model to unify the benchmark, and establish the absolute water level elevation of the surface water body.

[0122] like Figure 11 As shown, the hardware implementation of this system mainly relies on the GNSS-R water level monitoring terminal. This terminal integrates a dual-polarized antenna, radio frequency front-end, and digital receiver to support the operation of the sensing end signal acquisition module. The built-in edge computing unit is equipped with a water level inversion algorithm, which is used to execute the operation logic of the mirror point spatial positioning module, signal feature calculation module, and vertical ranging inversion module. At the same time, the local storage unit equipped with the terminal stores the station reference parameters for the absolute water level calculation module to call.

[0123] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for real-time surface water level inversion based on GNSS-R, characterized in that, Includes the following steps: S1: Detect direct and reflected signals from L-band GNSS satellites, simultaneously collect signal-to-noise ratio observations, carrier phase observations, and satellite pseudorange data under multi-constellation joint observations, and generate a multi-band GNSS raw observation signal set; S2: Based on the analysis of the transmission time of the original observation signal set of the multi-band GNSS, the orbital parameters are called to calculate the geometric transmission path, and the theoretical reflection position is screened using the digital elevation model and water body mask to obtain the coordinates of the water surface mirror reflection point; S3: For the region pointed to by the coordinates of the water surface mirror reflection point, calculate the correlation value between the reflection component and the copy code of the multi-band GNSS raw observation signal set, perform coherent integration and extract the peak delay, and obtain the delay Doppler map features of the reflection signal; S4: Call the delayed Doppler image feature of the reflected signal to read the time difference value, combine the coordinates of the water surface mirror reflection point to calculate the path difference and elevation angle sine value, use atmospheric, tidal and wind wave parameters to compensate for the foundation height, and establish the vertical distance value from the antenna phase center to the water surface. S5: Obtain the geodetic elevation parameters of the antenna reference ellipsoid based on the vertical distance from the antenna phase center to the water surface, calculate the difference between the geodetic elevation parameters of the reference ellipsoid and the vertical distance value, and generate the absolute water level elevation of the surface water body.

2. The real-time surface water level inversion method based on GNSS-R according to claim 1, characterized in that, The multi-band GNSS raw observation signal set includes signal-to-noise ratio observation sequence, carrier phase observation sequence, and multi-constellation pseudorange observation data; The coordinates of the water surface mirror reflection point include the geodetic coordinates of the mirror reflection point, the validity marker of the reflection point, and the boundary data of the water body mask. The delayed Doppler graph features of the reflected signal include the time delay at the specular reflection point, the Doppler frequency shift at the specular reflection point, and the peak power of the delayed Doppler graph; The vertical distance from the antenna phase center to the water surface includes the two-way propagation path difference, atmospheric and geophysical corrections, and antenna height inversion value. The absolute water level elevation of the surface water body includes the station ellipsoid height, the absolute water surface elevation, and the water level time series.

3. The real-time surface water level inversion method based on GNSS-R according to claim 1, characterized in that, The satellite pseudorange data is a distance metric obtained by multiplying the signal propagation time by the speed of light, and is used for timing alignment of the signal reception time. Digital elevation models are data containing surface elevation information, used to perform elevation anomaly correction on theoretical specular reflection points; The water mask is used to distinguish between land and water geographical data, determine the reflection point landing area attributes, and remove invalid areas.

4. The real-time surface water level inversion method based on GNSS-R according to claim 1, characterized in that, The peak delay is the maximum power point extracted from the two-dimensional delayed Doppler spectrum, used to characterize the time interval between the arrival of the reflected and direct signals; The vertical distance value is the vertical height from the antenna phase center to the water surface, obtained after time interval conversion, geometric path correction, and environmental parameter compensation. The reference ellipsoid geodetic elevation parameters are obtained by calling the receiver's precise single-point positioning results and filtering and smoothing them to determine the absolute height of the station antenna phase center relative to the reference ellipsoid surface.

5. The method for real-time surface water level inversion based on GNSS-R according to claim 1, characterized in that, The steps for obtaining S1 are as follows: The system detects right-hand circularly polarized direct signals emitted by L-band GNSS satellites and left-hand circularly polarized reflected signals reflected by the water surface, separates electromagnetic wave energy in different polarization directions, locks the signal frequency, and obtains the basic signal flow of the dual-polarization channel. The carrier and ranging code information of the basic signal stream of the dual-polarized channel are called to demodulate the signal-to-noise ratio (SNR) and carrier phase observations of the direct and reflected links. The validity of the SNR and carrier phase observations is screened to obtain the observation characteristic parameters of the direct and reflected links. Based on the observation characteristic parameters of the direct and reflected links, satellite pseudorange data under multi-constellation joint observation are collected synchronously. The satellite pseudorange data and characteristic parameters are time-aligned according to the reception time to establish a multi-band GNSS raw observation signal set.

6. The method for real-time surface water level inversion based on GNSS-R according to claim 1, characterized in that, The steps for obtaining S2 are as follows: Based on the multi-band GNSS raw observation signal set, the satellite signal transmission time is analyzed, the precise orbital parameters and satellite clock error parameters are called to correct the satellite clock deviation, the three-dimensional orbital position of the satellite's center of mass at the transmission time is calculated, and the instantaneous spatial position parameters of the satellite are obtained. The instantaneous spatial position parameters of the satellite and the preset geodetic coordinate parameters of the receiver are used to construct a signal transmission geometric link. The theoretical position where the signal propagation path is tangent to the reference ellipsoid is calculated based on the principle of geometric optics reflection, and the spatial vector of the theoretical specular reflection point is generated. The spatial vector of the theoretical mirror reflection point is used to introduce digital elevation model data for elevation anomaly correction. Water body mask data is used to determine the surface attributes of the reflection point's landing area and remove invalid areas to establish the coordinates of the water surface mirror reflection point.

7. The method for real-time surface water level inversion based on GNSS-R according to claim 1, characterized in that, The steps for obtaining S3 are as follows: For the area pointed to by the coordinates of the water surface mirror reflection point, the reflected signal data segment is extracted from the multi-band GNSS raw observation signal set, a local signal copy with the same pseudo-random code structure as the satellite transmission signal is generated, and the reflected signal components and the local copy code sequence are obtained. The reflected signal component and the local copy code sequence are called to construct a two-dimensional search grid in the time delay domain and Doppler frequency domain. Cross-correlation calculation is performed on the reflected signal component and the local copy code sequence and coherence integral is accumulated to obtain the signal power response intensity of each grid node and generate a two-dimensional delayed Doppler power spectrum array. Based on the two-dimensional delayed Doppler power spectrum array, the power values ​​of each grid node are traversed, the index of the location of the point with the maximum power is retrieved, and the time delay parameter is extracted as a feature quantity to obtain the delayed Doppler map feature of the reflected signal.

8. The real-time surface water level inversion method based on GNSS-R according to claim 1, characterized in that, The steps for obtaining S4 are as follows: The delayed Doppler map features of the reflected signal are called up, the time difference between the direct signal and the reflected signal arriving at the receiver is read, and a geometric link model is constructed by combining the coordinates of the water surface mirror reflection point and the satellite spatial coordinates. The geometric difference of the two-way path and the sine of the incident elevation angle of the satellite relative to the reflection plane are calculated to obtain the dual-station geometric delay and angle dataset. The time difference is converted into a distance by calling the bi-station geometric delay and angle dataset and introducing the vacuum speed of light constant. The geometric height is calculated by using the bi-path geometric difference and the sine of the incident elevation angle. The geometric path is corrected by combining the Earth curvature parameter to generate the geometric inversion foundation vertical height. Based on the geometric inversion foundation vertical height, atmospheric delay environment calibration parameters, tidal environment influence parameters, and wind and wave environment influence parameters are obtained. The error superposition and cancellation value of the foundation vertical height is calculated and dynamically compensated, and the vertical distance value from the antenna phase center to the water surface is established.

9. The method for real-time surface water level inversion based on GNSS-R according to claim 1, characterized in that, The steps for obtaining S5 are as follows: The observation epoch is determined based on the vertical distance from the antenna phase center to the water surface. The sequence of precise single-point positioning results from the receiver is retrieved. Points with jumps in positioning coordinates are removed by filtering and smoothing, and the average vertical height is calculated to obtain the reference geodetic elevation parameters of the station. The difference between the station's reference geodetic elevation parameters and the vertical distance from the antenna phase center to the water surface is calculated by calling the station's reference geodetic elevation parameters and the vertical distance from the antenna phase center to the water surface. The attitude error of the difference is corrected by the equipment installation tilt angle and signal quality factor, and a weighted calculation is performed to generate a dynamic elevation settlement value for the water surface. Based on the calculated dynamic elevation values ​​of the water surface, the parameters of the geoid refinement model for this region are introduced to calculate the elevation anomaly separation between the reference ellipsoid and the geoid. The water surface elevation is then converted to an orthographic system and unified with the benchmark to establish the absolute water level elevation of the surface water body.

10. A real-time surface water level inversion system based on GNSS-R, characterized in that, The system is used to execute the GNSS-R-based real-time surface water level inversion method according to any one of claims 1-9, comprising: a sensing end signal acquisition module, a mirror point spatial positioning module, a signal feature calculation module, a vertical ranging inversion module, and an absolute water level calculation module; The sensing end signal acquisition module is used to capture L-band signals through a dual-polarized antenna, acquire signal-to-noise ratio and carrier phase data through a radio frequency front-end and a digital receiver, synchronize multiple constellation pseudoranges, and establish a multi-band GNSS raw observation signal set. The mirror point spatial positioning module is used to call the original observation signal set of the multi-band GNSS to analyze the transmission time, calculate the geometric path in combination with the precise orbit parameters, use the digital elevation model to screen the landing area, and generate the coordinates of the water surface mirror reflection point. The signal feature calculation module is used to drive the edge computing unit to call the multi-band GNSS raw observation signal set to extract data for the area pointed to by the coordinates of the water surface mirror reflection point, generate a local copy code, perform cross-correlation calculation in the time delay and Doppler domain and accumulate coherent integrals to construct a two-dimensional power spectrum, traverse the grid nodes to retrieve the point with the maximum power, extract the time delay parameters, and obtain the delayed Doppler map features of the reflected signal. The vertical ranging inversion module is used to read the time difference by utilizing the delayed Doppler image features of the reflected signal, calculate the geometric path difference by combining the coordinates of the water surface mirror reflection point, introduce environmental parameters to compensate for the foundation height, and generate the vertical distance value from the antenna phase center to the water surface. The absolute water level calculation module is used to retrieve locally stored station reference parameters based on the vertical distance from the antenna phase center to the water surface, calculate the elevation difference, introduce a unified benchmark from the geoid model, and establish the absolute water level elevation of the surface water body.

Images