River water level inversion method and system considering altimeter footprint heterogeneity

Abstract

The application discloses a river water level inversion method considering altimeter footprint heterogeneity. The method comprises the following steps: acquiring a SAR altimeter observation waveform; calculating spatial geometric parameters through footprint gridding and observation geometry modeling; constructing a water body classification scene based on a water body mask and a river network, and classifying water bodies into main river channels, floodplains and adjacent isolated water bodies; performing electromagnetic scattering numerical integration based on differentiated surface roughness parameters to generate simulated echo waveforms; realizing waveform re-tracking and parameter inversion through coarse grid search and local optimization solution; calculating absolute water level and performing geophysical correction; performing multi-index quality control based on the interpretation rate of the simulated waveforms and the residual whiteness; and finally outputting water level products. The application overcomes the defects of the prior art, such as insufficient inversion precision and robustness in a complex surface mixed scene, and significantly improves the accuracy, reliability and physical interpretability of inland river water level monitoring.

Description

Technical Field

[0001] This invention relates to the field of satellite radar altimetry, and more specifically, to a river level inversion method that takes into account the heterogeneity of altimeter footprints and is applicable to complex inland water environments. Background Technology

[0002] River water level is a core parameter characterizing hydrological dynamics such as river flow volume and velocity, and it plays an irreplaceable role in flood forecasting, drought monitoring, water resource management and water conservancy project scheduling, as well as hydrological response research to global climate change. Traditional water level monitoring mainly relies on surface hydrological station networks. Although these networks can provide high-precision, high-frequency water level data, their spatial representativeness is limited, making it difficult to cover the vast global river systems, especially in remote mountainous areas, transboundary rivers, ecologically sensitive areas, and underdeveloped regions. The construction and long-term maintenance costs of these station networks are also high, and the number of global hydrological stations has been declining in recent years, further highlighting the inadequacy of regional and large-scale hydrological monitoring capabilities.

[0003] Against this backdrop, satellite radar altimetry, with its global coverage, periodic observations, and high consistency, has become an important means of supplementing and extending ground-based monitoring networks. This technology is based on the principle of active radar altimetry: a satellite platform transmits radar pulses towards the Earth's surface, receives the reflected echoes, and accurately determines the round-trip time from the satellite to the water surface by retracking the echo waveforms. Combined with precise orbit and attitude data, as well as a series of geophysical corrections for atmospheric, ionospheric, and solid tide data, the absolute elevation of the water surface, i.e., the water level, is ultimately obtained through inversion. Since Brown proposed the classic ocean echo model in the 1970s, radar altimetry has matured and been widely applied in sea level measurement. In 1998, Raney introduced delay-Doppler (DDoL) processing into the field of altimetry, significantly improving the track-direction resolution and signal-to-noise ratio of altimetry through along-track synthetic aperture technology. Subsequently, satellite platforms such as Sentinel-3 / SRAL significantly improved the resolution along the orbit through synthetic aperture processing, creating an opportunity to extend altimetry technology to inland water bodies with complex terrain and fragmented targets.

[0004] However, existing satellite radar altimetry systems were primarily designed for homogeneous, open marine environments. When applied to monitoring inland rivers with complex topography and diverse surface cover types, they face fundamental challenges due to the heterogeneity of observation geometry and surface scattering characteristics. The instantaneous illumination footprint of a radar altimeter is large in the trans-orbital direction (often on the order of kilometers), while the width of the target river is often only tens to hundreds of meters, much smaller than the footprint range. Therefore, within an observation footprint, the radar beam inevitably covers the main channel of the target river, the vegetation zones on both banks, seasonally flooded mudflats and wetlands, and adjacent still water bodies (such as lakes, ponds, and reservoirs), among other surface types. These different types of land features exhibit distinctly different electromagnetic scattering mechanisms under near-normal incident radar wave illumination: open, calm water surfaces primarily exhibit strong specular reflection; narrow, flowing river sections exhibit quasi-spectral scattering after attenuation by disturbed water surfaces; floodplains and mudflats exhibit diffuse scattering from rough surfaces; and vegetation exhibits a mixture of volume scattering and surface scattering. These scattered echoes, with their different physical mechanisms, superimpose upon reception, resulting in highly variable waveforms that often exhibit complex characteristics such as multiple peaks, broadening, and distorted leading edges. This deviates significantly from the classic Brownian model based on the assumption of a uniform sea surface. In this situation, each peak in the waveform often corresponds to a different physical scattering source. Accurately identifying and separating the effective signal component originating solely from the main channel of the target river from these overlapping peaks becomes the core challenge in inland water level inversion.

[0005] For water level extraction from complex waveforms, existing methods mainly develop along two technical paths:

[0006] (1) Empirical re-tracking methods based on waveform statistical features, such as OCOG (centroid shift method) and thresholding methods. These methods directly extract mathematical features from the waveform, which is simple to calculate and efficient. However, they mainly rely on the overall statistical characteristics of the waveform and lack consideration of the physical mechanism of waveform formation. In scenarios with significant footprint heterogeneity and severe waveform distortion, their stability and inversion accuracy drop sharply, and it is difficult to explain the physical causes of phenomena such as waveform multi-peaks and deformation.

[0007] (2) Physical retracking methods based on parametric models, such as the improved Brown model and the SAMOSA model specifically developed for synthetic aperture radar altimeters. These methods fit waveforms by constructing parametric mathematical models, and their adaptability is better than purely empirical methods. However, they are essentially still within the scope of mathematical curve fitting, and there is a lack of clear and universal physical correspondence between model parameters and surface physical properties. More importantly, existing models usually treat the entire footprint as a homogeneous or simply layered scatterer, failing to explicitly incorporate the spatial distribution of land cover types within the footprint and their differentiated scattering mechanisms. Therefore, when the waveform has multiple peaks, existing methods cannot reliably determine which peak corresponds to the target water body, and in practice, the strongest or sharpest peak is often simply identified as the water body peak. However, numerous observations show that nearby calm water bodies or wet mudflats often produce stronger mirror or mirror-like echoes than narrow, turbulent main channels, resulting in significant systematic biases or even complete errors in water level inversion.

[0008] In summary, the mixing of landform types within the radar altimeter's illumination footprint and the resulting heterogeneity in scattering characteristics are the most critical bottleneck currently hindering the improvement of accuracy and reliability of synthetic aperture radar (SAR) altimeter technology for river water level monitoring. To overcome this bottleneck, a novel technological approach is urgently needed to enhance the accuracy, robustness, and physical interpretability of satellite radar altimeter water level retrieval in complex inland scenarios. Summary of the Invention

[0009] This invention aims to overcome the shortcomings of existing technologies and provide a river water level inversion method that takes into account the heterogeneity of altimeter footprints. This method integrates multi-source spatial prior data to construct a refined water body scenario and establishes differentiated numerical echo models based on electromagnetic scattering mechanisms, thereby achieving high-precision and high-reliability river water level inversion from complex multi-peak waveforms through physical-driven methods.

[0010] To achieve the above objectives, the present invention adopts the following technical solution:

[0011] Firstly, a method for river level inversion that takes into account the heterogeneity of altimeter footprints includes the following steps:

[0012] S1. Data Acquisition and Preprocessing: Acquire observation waveform data, satellite orbit parameters, attitude parameters, instrument parameters, geophysical correction data, and spatial prior data including river network, water body mask, and digital elevation model data from synthetic aperture radar altimeter; perform noise removal and normalization on the observation waveforms and map them to a unified time delay grid.

[0013] S2. Footprint Gridding and Observation Geometric Modeling: The radar ground illumination area corresponding to a single observation epoch is discretized into a grid; based on the satellite orbital parameters, attitude parameters, and the digital elevation model, the geometric parameters of each grid cell relative to the satellite are calculated, including slant range, incident angle, and azimuth angle;

[0014] S3. Water body scene construction: Within the footprint network determined in step S2, water body pixels are identified based on the water body mask; the water body pixels are classified according to their connectivity with the centerline of the river network and their spatial position on the digital elevation model. The classification results of the water body pixels are: main river channel water body, floodplain / tidal flat water body connected to the main river channel, and adjacent isolated water body not connected to the main river channel.

[0015] S4. Numerical echo simulation with electromagnetic scattering constraints: Different surface roughness parameters are set for the main channel water body and adjacent isolated water bodies respectively; Based on the geometric parameters calculated in step S2 and the classification results in step S3, electromagnetic scattering numerical integration is performed on the footprint grid to generate simulated echo waveforms; wherein, water level changes are simulated by changing the elevation of the grid cells corresponding to the main channel water body, and this elevation change is converted into the overall translation of the simulated waveform on the time delay axis through geometric relationships;

[0016] S5. Waveform re-tracking and parameter inversion: Construct an objective function that minimizes the difference between the simulated waveform and the observed waveform; First, perform a coarse grid search in the interval near the initial time delay value determined based on the digital elevation model to obtain the initial value of the inversion parameters; Then, perform local optimization solution starting from the initial value to obtain the optimal water level time delay parameters;

[0017] S6. Distance and Elevation Calculation and Standard Correction: Based on the optimal water level time delay parameters, calculate the slant distance and initial ellipsoidal height, and apply wet / dry tropospheric correction, ionospheric correction, and solid tide correction one by one, and convert to the absolute water level under the unified geoid reference;

[0018] S7. Multi-index quality control: Calculate the explanatory power of the simulated waveform to the observed waveform, and the whitening degree of the fitting residual; perform quality judgment and screening on the inversion results based on the explanatory power and the whitening degree;

[0019] S8. Product Output: Output a water level product that includes the absolute water level, inversion uncertainty, and quality indicator.

[0020] Furthermore, in step S3, the classification of water body pixels is specifically as follows: water body pixels that are connected to the center line of the river network and located within the main channel space are identified as main river channel water bodies; water body pixels that are connected to the main river channel water bodies but located outside the main channel space are identified as floodplain or tidal flat water bodies; and water body pixels that exist within the footprint grid but have no topological connection with the main river channel water bodies are identified as adjacent isolated water bodies.

[0021] Furthermore, in step S4, the surface roughness parameter is a mean square slope parameter based on a geometrical optical approximation model, used to characterize the statistical properties of the microscopic undulations of the water surface; wherein, the mean square slope parameter value set for adjacent isolated water bodies is less than the mean square slope parameter value set for the main river channel water body.

[0022] Furthermore, the electromagnetic scattering numerical integration in step S4 specifically includes: for each water body grid cell, calculating its backscattering coefficient according to its category by calling the corresponding surface roughness parameter, and calculating the range attenuation factor and antenna pattern gain weight in the time delay-Doppler domain by combining its geometric parameters, and generating the component waveform of this type of water body through numerical integration; and performing incoherent linear superposition of the component waveforms of the main channel water body, the adjacent isolated water body, and a weak background noise term to generate the simulated echo waveform.

[0023] Furthermore, in step S5, the range of the coarse grid search is determined based on the prior water surface elevation provided by the digital elevation model, combined with the time delay reference value obtained by converting the satellite orbital altitude, and fluctuating up and down by a preset range of physical water level changes, wherein the range of physical water level changes is 2 meters to 10 meters.

[0024] Furthermore, in step S5, the nonlinear local optimization solution adopts the Levenberg-Marquardt algorithm, and during the optimization process, a value range constraint based on physical experience is applied to the surface roughness parameter to prevent the parameter optimization from falling into non-physical solutions.

[0025] Furthermore, step S5 also includes: calculating and outputting the inversion uncertainty corresponding to the water level delay parameter based on the posterior covariance information obtained from the optimization process.

[0026] Furthermore, in step S7, the quality determination is as follows: when the explanatory rate is not lower than the first threshold and the whitening degree is not lower than the second threshold, the epoch inversion result is determined to be of qualified quality; otherwise, it is discarded or marked. Wherein, the first threshold is 0.75 and the second threshold is 0.6.

[0027] Furthermore, step S8 also includes: fusing multiple qualified water level products obtained by inverting the same virtual station from different satellite orbits, and using orbital intersection consistency constraints or ground reference water level data to perform system deviation correction, thereby generating a time-series continuous river water level product.

[0028] Secondly, a river level inversion system that takes into account the heterogeneity of altimeter footprints includes:

[0029] The data access and preprocessing module is used to acquire and preprocess the observation waveform data, satellite orbit parameters, attitude parameters, instrument parameters, geophysical correction data and space prior data of the synthetic aperture radar altimeter;

[0030] The geometric modeling module is used to discretize the radar ground illumination area into a grid and calculate the geometric parameters of each grid cell.

[0031] The scene construction module is used to identify and classify water body pixels within the footprint grid, mainly river water bodies, floodplain / tidal flat water bodies, and adjacent isolated water bodies.

[0032] The numerical echo simulation module is used to perform numerical integration of electromagnetic scattering on the footprint grid based on differentiated surface roughness parameters and geometric parameters, and generate simulated echo waveforms.

[0033] The parameter inversion module is used to solve for the optimal water level delay parameters that minimize the difference between the simulated waveform and the observed waveform through coarse grid search and local optimization.

[0034] The solution and correction module is used to calculate the absolute water level based on the optimal parameters and perform geophysical corrections.

[0035] The quality control module is used to judge and screen the inversion results based on the interpretation rate of the simulated waveform and the whitening degree of the fitting residuals.

[0036] The product output module is used to output water level products that include absolute water level, inversion uncertainty, and quality indicators.

[0037] The advantages of this invention over the prior art are:

[0038] This invention significantly improves the accuracy and robustness of water level inversion in complex inland water body scenarios: Addressing the challenges of multi-peak and deformable waveforms caused by the mixing of land cover types within the radar altimeter footprint, this invention constructs a refined water body scenario of "main channel - floodplain - adjacent still water body" by fusing multi-source spatial prior data, and establishes a differentiated numerical echo model based on electromagnetic scattering mechanisms. This model can physically drive the separation and interpretation of superimposed echo signals from different land cover types (such as disturbed narrow channels and adjacent calm ponds), enabling the inversion process to accurately identify and lock the effective signal components of the target main channel water surface. This significantly reduces the systematic bias caused by misidentifying the strongest peak as a water body peak, resulting in a substantial improvement in the accuracy and reliability of water level inversion in complex environments.

[0039] This invention provides explicit physical interpretability, enhancing the transparency and credibility of the inversion process: Compared to traditional empirical re-tracking algorithms or parametric fitting models, the "numerical integration-component synthesis" echo model constructed in this invention has a solid physical foundation in electromagnetic scattering. Model parameters (such as surface mean square slope and water level delay) directly correspond to surface physical properties, and water level changes drive the overall translation of the simulated waveform through geometric relationships. Simultaneously, physical constraints (such as roughness parameter ranges and DEM-based initial delay constraints) are introduced during the inversion process to prevent non-physical interpretations. This physical-driven modeling and inversion method not only improves the accuracy of the results but also makes the entire inversion process and results more interpretable and physically meaningful.

[0040] A systematic quality control and uncertainty quantification system was constructed to ensure product quality consistency. This invention innovatively proposes a joint quality control index system of "model explanation rate" and "residual whitening degree." Through dual-threshold gating and a weighted scoring mechanism, it can automatically identify and eliminate abnormal observation epochs affected by extreme interference (such as strong specular reflection or occlusion). Simultaneously, using posterior covariance information from the optimization process, the uncertainty of the water level inversion results is quantitatively calculated and output. This rigorous quality control, encompassing both goodness of fit and residual statistical characteristics, combined with quantitative assessment of uncertainty, ensures that the final water level products (including single-epoch results and fused time series) possess reliable quality standards and known confidence levels.

[0041] This invention achieves automated and scalable operational processing capabilities, improving monitoring efficiency and product usability. The method forms a complete and structured processing flow from data preprocessing, scenario construction, physical modeling, parameter inversion to product generation and quality control. This flow can be automated and adapted to rivers of varying widths (through adaptive grid resolution) and scenarios of varying complexity. Ultimately, through multi-track fusion and system bias correction, it can generate spatiotemporally continuous and benchmark-unified river water level time series products, directly serving operational applications and model assimilation in hydrological monitoring, flood forecasting, and water resource management, significantly expanding the practical value and application scope of satellite radar altimetry in inland water monitoring. Attached Figure Description

[0042] To more clearly illustrate the technical solutions in the embodiments of this application, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0043] Figure 1 This is a flowchart of the overall process of a river water level inversion method that takes into account the heterogeneity of altimeter footprints according to the present invention.

[0044] Figure 2 This is a simplified schematic diagram of the observation geometry in an embodiment of the present invention.

[0045] Figure 3 This is a schematic diagram of the water scene construction according to an embodiment of the present invention.

[0046] Figure 4 This is a flowchart of the numerical echo simulation according to an embodiment of the present invention.

[0047] Figure 5 This is a flowchart of waveform re-tracking and parameter inversion according to an embodiment of the present invention.

[0048] Figure 6 This is a flowchart of distance and elevation calculation and geophysical correction according to an embodiment of the present invention.

[0049] Figure 7 This is a multi-index quality control judgment logic diagram according to an embodiment of the present invention.

[0050] Figure 8 This is a flowchart of the product output and time series fusion process according to an embodiment of the present invention.

[0051] Figure 9 This is a comparison diagram of numerical simulation waveforms and satellite measured waveforms in a typical river section according to an embodiment of the present invention.

[0052] Figure 10This is a time series comparison chart of the water level inversion results of the method of this invention and the existing mainstream retracking algorithms.

[0053] Figure 11 These are numerical simulation waveforms for four typical multi-water-body interference scenarios.

[0054] Figure Labels: 100 Data Access and Preprocessing Module, 101 Observation Data Input, 102 Satellite Orbit Data, 103 Satellite Attitude Data, 104 Instrument Parameters, 105 Geophysical Correction Data, 106 Space Priors, 107 Noise Basis Removal, 108 Unified Delay Grid, 200 Geometric Modeling Module, 201 Origin (O), 202 Track Axis (x-axis), 203 Lateral Track Axis (y-axis), 204 Vertical Axis (z-axis), 205 Satellite Platform (S), 206 Satellite Velocity Vector (v), 207 Footprint Grid Element (F), 208 Satellite-Ground Slant Range Vector (R), 209 Angle of Incidence (θ), 210 Azimuth (φ), 221 Satellite Altitude (Hs), 300 Scene Construction Module, 301 Footprint Grid, 302 River Centerline, 303 Main Channel Water, 304 Adjacent Water (Not Connected to Main Channel), 305 Floodplain / Tidal Flat Water, 306 Water Mask, 400 Numerical Echo Simulation Module, 402 Mapping Matrix, 403 System Point Target Response (PTR), 404 Antenna Pattern Gain, 405 Fresnel Reflection Coefficient, 406 Surface Slope Statistics, 407 Water Scattering Model, 408 Numerical Integrating Unit, 409 Doppler Appearance Aggregation, 410 Simulated Waveform, 411 Background Noise Term, 412 Component Waveform Generation, 413 Window Function, 414 Sampling and Quantization, 415 Linear Superimposed Unit, 416 Main Channel Scattering Model, 417 Floodplain Scattering Model, 418 Adjacent Water Scattering Model, 419 420 Main channel integration, 421 Floodplain integration, 422 Adjacent water volume integration, 423 Main channel component waveform, 424 Floodplain component waveform, 425 Adjacent water volume component waveform, 500 Parameter inversion module, 501 Objective function, 502 Coarse grid search, 503 Initial parameters, 504 LM local optimization, 505 Physical constraints / regularization, 506 Convergence determination, 507 Optimal parameters, 508 Jacobian calculation, 509 Uncertainty assessment, 510 Observed waveforms, 511 Numerical model interface, 600 Solution and correction module, 601 Round trip delay, 602 Slant distance calculation, 603 Water surface ellipsoid elevation (uncorrected), 604 Wet troposphere correction, 605 Dry troposphere correction, 606 Ionospheric correction. 607 Fixed Tide Correction, 608 Load Tide Correction, 609 Polar Tide Correction, 610 Dynamic Atmospheric Correction, 611 Correction Summary, 612 Corrected Ellipsoid Elevation, 613 Geoid Model, 614 Unified Baseline Water Level, 615 Satellite Ellipsoid Elevation, 700 Quality Control Module, 701 Model Interpretation Rate.702 Residual Whitening, 703 Gated Merging Node, 704 Weighted Scorer, 705 Comprehensive Judgment, 706 Reconstruction Trigger / Failure, 707 Quality Control Pass, 708 Explanation Rate Threshold Gating, 709 Whitening Threshold Gating, 800 Product Output Module, 801 Virtual Station Water Level, 802 Water Level Time Series, 803 Epitome Uncertainty, 804 Quality Marker, 805 Multitrack Fusion, 806 Offset / Scale Correction, 807 Reference Water Gauge / Control Point, 808 Crosspoint Constraint, 809 Temporal Consistency Screening / Outlier Removal, 810 Product Archiving and Release, 811 Metadata and Version Management, 812 Time Resampling / Interpolation, Detailed Implementation

[0055] In the following description, specific details such as particular system architectures and techniques are set forth for illustrative purposes and not for limiting purposes, in order to provide a thorough understanding of the embodiments of this application. However, those skilled in the art will understand that this application may also be implemented in other embodiments without these specific details.

[0056] Reference Figure 1 , Figure 1 A schematic diagram of the overall process of the inland water level physical inversion method for the Sentinel-3 SAR model according to an embodiment of the present invention is shown. The method mainly includes data access and preprocessing, footprint gridding and observation geometry modeling, water body scene construction, numerical echo simulation with electromagnetic scattering constraints, waveform retracking and parameter inversion, distance and elevation calculation and standard correction, multi-index quality control and product output, etc.

[0057] The specific implementation of the river water level inversion method of the present invention, which takes into account the heterogeneity of altimeter footprints, includes the following core steps, and its detailed technical features are described below:

[0058] S1. Data Acquisition and Preprocessing

[0059] In the data access and preprocessing stage, observational data from the Sentinel-3 satellite (mainly 20 Hz waveform data) is input 101, and the accompanying satellite orbit data 102 and satellite attitude data 103 are read simultaneously. At the same time, instrument parameters 104 are loaded, specifically including antenna pattern gain 404, system point target response 403, pulse repetition frequency, and echo sampling rate, among other system configurations. Furthermore, geophysical correction data 105 and spatial priors including GRWL river networks, SWOT water surface elevation, or DEM are introduced 106. After noise floor removal 107, the raw data is aligned in the temporal / spatial domain, mapped to a unified time-delay raster 108, and a standard data cube required for inversion is established.

[0060] S2. Footprint Gridding and Observational Geometric Modeling

[0061] This step constructs an accurate observation geometry model for each 20 Hz observation epoch. First, the ground area (footprint envelope) illuminated by the radar beam at that epoch is discretized into a high-resolution grid. The grid resolution is adaptively adjusted according to the target river width; for example, a resolution of 5 meters is used for narrow channels (<100 meters), 15 meters for medium channels (100-500 meters), and 30 meters for wide channels (>500 meters). Based on the satellite's real-time orbital position (longitude, latitude, ellipsoidal height) and attitude angle, combined with the WGS84 Earth ellipsoid model and the digital elevation model (DEM) loaded in step S1, the spatial geometric relationship between the center point of each grid cell and the phase center of the satellite antenna is calculated.

[0062] Specifically, for each grid cell, the following parameters are calculated: Slant range R: The straight-line distance from the satellite to the grid cell. Incident angle θ: The angle between the slant range direction and the grid cell normal direction (calculated from the DEM). Azimuth angle φ: The angle between the projection of the slant range vector onto the horizontal plane and the satellite's trajectory direction (flight direction).

[0063] like Figure 2 As shown, a local rectangular coordinate system is established to demonstrate the spatial geometric relationship between the satellite observation platform and the Earth's surface footprint grid. The geometric model uses the Earth's surface reference point as the origin 201 (O) and constructs a three-dimensional orthogonal coordinate system, including the track axis 202 (x-axis) indicating the direction of the satellite's flight trajectory, the transverse track axis 203 (y-axis) perpendicular to the flight direction, and the vertical axis 204 (z-axis) pointing towards the zenith.

[0064] Satellite platform 205(S) is located on the vertical axis, and its vertical distance relative to the origin 201 is represented by satellite altitude 221(Hs). The satellite moves in the direction indicated by satellite velocity vector 206(v). On the surface plane, footprint grid element 207(F) represents a discrete ground scattering unit. The line connecting satellite platform 205 and footprint grid element 207 forms satellite-to-ground slant range vector 208(R), which is the geometric basis for radar echo delay calculation.

[0065] Based on the above geometric elements, two key angular parameters that determine the electromagnetic scattering characteristics are defined: (1) Incident angle 209 (θ): defined as the angle between the vertical axis 204 and the star-to-ground slant distance vector 208. This angle directly determines the specular reflection intensity in the calculation of the backscattering coefficient. (2) Azimuth angle 210 (φ): defined as the angle between the trace axis 202 and the star-to-ground slant distance vector projected onto the ground surface in the xOy plane. This angle is used to determine the calculation direction of the Doppler frequency shift.

[0066] Simultaneously, based on the orientation of the grid cell relative to the antenna, the antenna pattern gain function is queried to obtain the gain weight G(θ, φ) of that cell. This step constructs an accurate mapping matrix from geographic space (grid cell) to observation space (time-delay-Doppler cell), which forms the basis for subsequent numerical integration.

[0067] S3. Water Body Scene Construction

[0068] Within the footprint network determined in step S2, water body pixels are identified based on the water body mask; the water body pixels are classified according to their connectivity with the centerline of the river network and their spatial position on the digital elevation model. The classification results of the water body pixels are: main river channel water bodies, floodplain / tidal flat water bodies connected to the main river channel, and adjacent isolated water bodies not connected to the main river channel.

[0069] This invention establishes only the water scattering scene within the radar footprint envelope at each 20 Hz observation epoch, without separately modeling non-water features. The physical basis lies in the combined constraints of time delay-Doppler processing and antenna pattern and time delay window, concentrating effective echo energy within the main lobe and target time delay range. Near-normal incident water surfaces scatter primarily in a near / quasi-mirror manner, constituting the dominant echo component, while vegetation, moist soil, and bare land primarily scatter diffusely, and their energy contribution to the same main peak gate is usually negligible. Specifically, the GRWL river network and centerline, GSWO water frequency, and elevation data (preferably using SWOT water surface elevation, defaulting to SRTM DEM data) are overlaid within the footprint network. A binary water mask is constructed with GSWO occurrence frequencies ≥10%. After morphological opening and closing, connected component filtering (removing isolated patches smaller than 3–5 pixels), and hole filling processing, it is cropped to the footprint range and resampled to grid resolution. Based on their connectivity with the GRWL centerline and spatial location, water bodies within the footprint are divided into three categories: (i) those connected to the centerline and located within the main channel zone are defined as the main channel; (ii) those located within the footprint but not within the main channel zone are defined as floodplains / tidal flats; and (iii) those not connected to the main channel are defined as adjacent water bodies (lakes / reservoirs / ponds). Non-water body pixels do not participate in subsequent numerical integration calculations. This processing creates a fractionalized scene within the footprint domain of "main channel—floodplain / tidal flat—adjacent still water body," with clear boundaries, physical consistency, and numerical stability.

[0070] like Figure 3 , Figure 3 The results of surface classification and its spatial topological relationships within the radar observation field of view are presented. Figure 3The dashed elliptical region represents the effective footprint network 301 of the radar altimeter, which defines the geographical boundary of the numerical integration. Within the footprint network 301, based on multi-source spatial prior data, binarized water masks 306 (corresponding to all marked non-blank areas in the figure) were identified and extracted, and then subjected to refined object-level classification.

[0071] (1) River centerline 302: extracted from prior databases such as GRWL, indicating the geometric skeleton direction of the target water flow.

[0072] (2) Main channel 303: A strip of continuous water distributed along the center line of the river, which is the main target of the water level inversion in this invention.

[0073] (3) Adjacent water bodies 304: Independent still water bodies (such as oxbow lakes, ponds, reservoirs, etc.) that are distributed on both sides of the main channel and are not spatially connected to the main channel. These water bodies usually appear as high-brightness, strong scattering sources and are the main sources of interference that produce false peaks.

[0074] (4) Floodplain / Tidal Flat 305: The edge area that has spatial connectivity with the main channel but may be exposed or covered by vegetation during the dry season. The scattering characteristics of this area vary with the rise and fall of the water level and usually exhibit mixed scattering characteristics.

[0075] Based on a unified time-delay grid 108 and spatial priors 106, refined land cover classification is performed within the satellite footprint network 301. Specifically, the main river channel 303, adjacent water bodies 304 not connected to the main river channel, and connected floodplains / tidal flats 305 are identified. A binary water body mask 306 containing each independent land cover object is generated accordingly, providing geometric boundary constraints for component modeling. Through this classification, the present invention deconstructs the complex inland surface into several independent units with different physical properties (such as connectivity and roughness), providing basic scene support for subsequent differentiated parameter modeling.

[0076] S4. Numerical echo model with electromagnetic scattering constraint

[0077] Different surface roughness parameters are set for the main channel water body and the adjacent isolated water body respectively; based on the geometric parameters calculated in step S2 and the classification results in step S3, electromagnetic scattering numerical integration is performed on the footprint grid to generate simulated echo waveforms; wherein, water level changes are simulated by changing the elevation of the grid cells corresponding to the main channel water body, and the elevation change is converted into the overall translation of the simulated waveform on the time delay axis through geometric relationships.

[0078] This step establishes a numerical echo model based on the radar electromagnetic scattering mechanism, employing a strategy of "infinite element scattering calculation—component waveform synthesis—full waveform incoherent superposition." The specific implementation process is as follows:

[0079] (1) For micro-elements within the grid that are identified as different water body types (such as main channels, floodplains, and adjacent still water ponds), a parameterized model of the backscattering coefficient is established based on the geometric optics approximation theory:

[0080]

[0081] The Fresnel reflection coefficient characterizes the dielectric properties of water, and the mean square slope (mss) is introduced as a statistical measure to characterize the microscopic roughness of the water surface. To realistically reflect the physical differences in echo amplitude and energy distribution among different water bodies in complex inland scenarios, this step sets differentiated roughness parameters for different independent water body components within the footprint: a higher initial mean square slope is set for the main river channel disturbed by water flow to simulate the attenuation of its quasi-specular scattering characteristics; a lower initial mean square slope is set for relatively closed and calm oxbow lakes or paddy fields to simulate their high-intensity echo characteristics approaching specular reflection.

[0082] (2) Couple the above physical scattering terms with the geometric and response characteristics of the satellite observation system to generate component power waveforms. Based on the real-time orbital state vector of the satellite and the three-dimensional spatial coordinates (longitude, latitude, and elevation) of each grid element, calculate the two-way transmission delay and Doppler frequency shift of each element relative to the phase center of the radar antenna. Perform weighted numerical integration on the footprint grid in the time-delay-Doppler domain, and apply multiple weighted modulations to the energy of each element, including the backscattering coefficient calculated by the geometric optics model, the inverse fourth-power geometric attenuation, the antenna pattern gain covering the main lobe and side lobe effects, and the on-track Doppler beam sharpening function of the SAR mode. Subsequently, treat the weighted element energy as an impact response and convolve it with the system point target response function of the radar altimeter (usually approximated as the Sinc square function) to independently generate the power waveforms of each water body component in the time-delay domain. Specific calculation formula:

[0083]

[0084] Where λ is the radar wavelength; R is the distance from the phase center of the satellite antenna to the pixel; G is the antenna pattern gain function, covering the main lobe and side lobe effects; W is the on-track Doppler beam sharpening function in SAR mode; and PTR is the system point target response function of the radar altimeter, which is approximately the Sinc squared function in this embodiment.

[0085] (3) Based on the incoherent characteristics of radar echoes, component waveform synthesis, sampling matching, and physical constraint construction are performed. The specific calculation formula is as follows:

[0086]

[0087] The total power of the synthesized waveform (Ptotal) is equal to the incoherent linear superposition of the main channel component waveform (Priver) and the adjacent interfering water body component waveform (Pint), and a weak background noise term with an energy ratio of no more than 5% is introduced to absorb sporadic residual non-water body scattering or thermal noise.

[0088] The synthesized total power waveform is resampled and truncated strictly according to the actual observation range gate resolution, generating a simulated echo consistent with the satellite measured data in terms of sampling structure. In this process, water level is explicitly defined as an independent variable that changes the grid elevation. The vertical change of water level is directly converted into the overall translation of echo time delay through geometric relationships, driving peak position movement and waveform morphology evolution. This mechanism establishes a deterministic physical constraint of "water level change - time delay translation - peak position movement", no longer relying on the empirical rule of "the strongest peak is the water body peak", thus reducing the systematic bias of the inversion.

[0089] like Figure 4 This diagram illustrates the complete signal processing chain and its logical relationships, from geophysical parameter input and component scattering modeling to the final radar waveform generation. The diagram uses a block flowchart to represent the forward modeling process based on physical mechanisms. This process calculates the scattering contributions of different water bodies in parallel, generates component waveforms through numerical integration, and simulates radar receiver processing after linear superposition. Specifically, it includes the following key processing nodes:

[0090] (1) Parameter input and mapping (404, 405, 406, 402):

[0091] 405 Fresnel Reflectance Coefficient and 406 Surface Slope Statistics: Input Parameters. These two parameters together determine the backscattering coefficient of a surface unit. Among them, the slope statistics 406 (such as root mean square slope) is used to characterize the roughness of the water surface and is the physical basis for distinguishing the scattering characteristics of still water and flowing water.

[0092] 402 Mapping Matrix: Establishes a spatial mapping relationship between the surface footprint grid and radar echo delay gates. It precisely quantifies each grid cell in geospatial space into an energy contribution weight for a specific distance gate.

[0093] (2) Differential scattering modeling (407, 408):

[0094] 407 Water Body Scattering Model: This model includes three independent physical scattering models: the main channel (416), the floodplain (417), and adjacent water bodies (418). It is designed for... Figure 3 Different surface types are identified and assigned different roughness parameters or scattering mechanisms (such as specular reflection or Bragg scattering) to accurately reflect the heterogeneity of real-world scenes.

[0095] 408 Numerical Integration Unit: Under the constraint of mapping matrix 402, spatial integration operations 419, 420, and 421 are performed on the three types of water bodies mentioned above. This step converts the spatially distributed scattered energy into a power sequence distributed along the time axis.

[0096] (3) Component waveform generation 412:

[0097] The output consists of three independent waveforms: the main channel (422), the floodplain (423), and the adjacent water body (424). These waveforms are not superimposed in the time domain and represent the original contribution of each water body to the total echo, providing direct observational evidence for analyzing interference sources (such as strong reflections from adjacent water bodies).

[0098] (4) Signal synthesis and processing (415, 409, 413, 414):

[0099] 415 Linear Superposition Unit: Based on the superposition principle of electromagnetic waves, the waveforms of each component are coherently or incoherently superimposed with the background noise 411. This step reassembles the decoupled physical model into the original signal of the entire scene.

[0100] 409 Doppler Appearance Aggregation: For delayed-Doppler (SAR) altimeter modes, beams from different viewpoints are weighted and aggregated to form multi-view waveforms.

[0101] Window function 413 and sampling quantization 414: Hardware characteristics of the analog radar receiver. Windowing smooths the signal sidelobes, and sampling quantization converts the analog signal into a digital signal, ultimately outputting an analog waveform 410 consistent with actual observations.

[0102] Forward physical modeling is performed based on the constructed scene. The main channel scattering model 416, the floodplain scattering model 417, and the adjacent water body scattering model 418 are invoked respectively. Combined with Fresnel reflection coefficient 405 and surface slope statistics 406 (such as root mean square slope), electromagnetic scattering integration is performed on the grid cells. The integration results of each component are incoherently superimposed using a linear superimposed unit 415, and then windowed / windowed processing 413 and sampling and quantization 414 are applied sequentially to generate a simulated waveform 410 consistent with the actual observed physical properties.

[0103] Through the above process, this invention achieves a physical mapping from "surface scene" to "radar signal". By utilizing the principle of linear superposition to deconstruct the complex scene into independently controllable component modules, not only is the physical fidelity of the simulation improved, but an interpretable forward model is also provided to support the subsequent inversion of the main channel water level through waveform decontamination or retracking algorithms.

[0104] S5. Waveform Re-tracking and Parameter Inversion

[0105] Construct an objective function that minimizes the difference between the simulated waveform and the observed waveform; first, perform a coarse grid search in the vicinity of the initial time delay value determined based on the digital elevation model to obtain the initial value of the inversion parameters; then, perform local optimization solution starting from the initial value to obtain the optimal water level time delay parameters.

[0106] Based on the physically driven simulated waveform generated by S4, this step employs a two-stage strategy of "coarse search - fine inversion." By minimizing the residual between the simulated waveform and the satellite-measured waveform, the accurate water level and related physical parameters are obtained through inversion. The specific implementation process is as follows:

[0107] (1) Constructing an objective function based on the least squares criterion

[0108] This step establishes an objective function based on the least squares criterion to quantify the degree of matching between the simulated and measured waveforms in the power domain. This objective function is defined as the sum of squared residuals between the simulated and measured waveforms at each sampling point, and includes regularization constraints to limit parameter deviations from the physically feasible range. The parameter vector to be inverted includes at least the zero-point delay offset, amplitude scaling factor, and surface roughness (mean square slope). The zero-point delay offset, as the primary control parameter, directly determines the overall translation of the simulated waveform on the time axis, corresponding to the precise correction of the radar wave's two-way transmission time. The amplitude scaling factor is used to correct the constant deviation in the backscattering coefficient model, and the surface roughness is used to adjust the peak sharpness and width characteristics of the waveform.

[0109] (2) Implement the first stage of coarse grid search to obtain robust initial values.

[0110] To address the common issues of non-convexity of the objective function and local minima in multi-peak waveform fitting, a coarse global search is first performed. Based on high-precision DEM elevation, historical hydrological statistics, or prior water levels from virtual stations, and combined with satellite orbital altitude conversion, a reference echo delay is obtained. This delay is used to define a physically reasonable search interval for zero-point offset (e.g., a delay range corresponding to a 10-meter fluctuation in DEM elevation). Simultaneously, a wide-range discrete grid is defined for nonlinear parameters such as amplitude scale and roughness. The objective function value is calculated by traversing this multi-dimensional parameter grid, searching for the parameter combination with the smallest residual, which is then used as the initial value for subsequent fine-tuning. This operation effectively prevents the algorithm from falling into erroneous convergence traps caused by strong non-target interference peaks, ensuring the global convergence of the inversion process.

[0111] (3) Implement the second stage of local optimization and perform fine estimation.

[0112] Within the initial value neighborhood determined by the coarse search, a nonlinear least squares optimization algorithm (such as the Levenberg-Marquardt algorithm) is used for iterative solution. During the iteration process, the partial derivatives of the objective function with respect to each parameter are calculated using numerical difference or analytical derivation to guide the parameters to converge quickly towards the optimal solution. To ensure the physical rationality of the results, physical boundary and regularization constraints are applied to the parameters during the optimization process: for the zero-point offset parameter of the time delay, a prior penalty term based on the DEM theoretical time delay is introduced to suppress jumps that deviate too much from the physical truth value; strict upper and lower limits are set for the surface roughness and amplitude parameters (such as limiting the roughness to positive values ​​and within the empirical range) to prevent non-physical "false solutions" or overfitting phenomena caused by noise.

[0113] (4) Evaluate uncertainty and output quality indicators

[0114] After iteratively converging to obtain the optimal zero-point offset for time delay, the river water level elevation (relative to a unified vertical reference) is converted into the real-time satellite orbital altitude and atmospheric propagation correction term. Simultaneously, based on the posterior covariance matrix of the parameters calculated by the inversion theory, the variance components corresponding to the time delay parameters are extracted and propagated to the water level domain. The standard deviation of the water level is calculated and output as an uncertainty index. The final output includes the accurately inverted water level value, inversion uncertainty, waveform fitting residuals, and algorithm convergence status indicators, providing complete quality information for subsequent hydrological applications or data assimilation.

[0115] like Figure 5 As shown, this demonstrates the complete iterative optimization process of using numerical models to approximate observational data in order to extract physical parameters.

[0116] This process employs a two-stage strategy of "coarse search + fine refinement." By constructing a differentiable objective function, it uses the Levenberg-Marquardt (LM) algorithm under physical constraints to find the optimal solution and finally evaluates the reliability of the results. Specifically, it includes the following key processing nodes:

[0117] (1) Input and interface construction (510, 511, 505):

[0118] Observed waveform 510: The echo signal actually received by the radar altimeter (i.e., Figure 4 The corresponding real data in the processing flow is the fitting target (Ground Truth) for parameter inversion.

[0119] Numerical Model Interface 511: Connection Figure 4 The numerical simulator's API is shown. During the inversion process, this API generates simulated waveforms in real time based on the parameter set of the current iteration, which are then compared with the observed waveforms.

[0120] Physical constraints / regularization 505: Parametric boundary conditions set based on prior knowledge (such as non-negative constraints on surface roughness). This module plays a regularization role in the optimization process, preventing the algorithm from converging to non-physical local minima under low signal-to-noise ratio or ill-conditioned conditions.

[0121] (2) Initialization and objective function (501, 502, 503):

[0122] Objective function 501: Defines a mathematical measure of the difference between the observed waveform and the simulated waveform (such as the maximum likelihood function or the least squares residual).

[0123] Coarse Grid Search 502: A global optimization strategy. It traverses a preset range within the parameter space with a large step size to calculate the objective function value.

[0124] Initial parameter 503: The optimal starting point determined by coarse-grid search. Providing high-quality initial values ​​is crucial to ensuring fast convergence of the nonlinear least squares problem and avoiding getting trapped in local optima.

[0125] (3) Iterative optimization and convergence (504, 506):

[0126] LM Levenberg-Marquardt Local Optimization (LM Local Optimization): The core iterative engine. It combines the advantages of gradient descent and Gauss-Newton methods, and continuously modifies parameters to minimize the objective function by dynamically adjusting the damping factor under the physical constraints.

[0127] Convergence Decision 506: Logical Decision Node (Diamond). Checks if the residual change or parameter update step size of the current iteration is less than a preset threshold. If the condition is met, convergence is determined and the result is output; otherwise, it returns to 504 to continue the next iteration.

[0128] (4) Result Calculation and Evaluation (507, 508, 509):

[0129] Optimal parameter 507: The final set of physical parameters output after iterative convergence, which usually includes the slant distance from the satellite to the water surface (corresponding to the water level), the root mean square slope of the water surface, the backscattering coefficient, etc.

[0130] Jacobi calculation 508: At the optimal parameter points, calculate the partial derivative matrix of the numerical model with respect to each inversion parameter. This matrix reflects the sensitivity of the model output to parameter variations.

[0131] Uncertainty assessment 509: Based on the Jacobian matrix and residual statistics, calculate the posterior covariance matrix of the parameters. This step quantitatively provides the confidence interval or standard deviation of the inversion results (such as water level), providing a reliable error measure for subsequent data assimilation or quality control.

[0132] A differentiable objective function 501 is constructed between the simulated waveform 410 and the observed waveform 510. The inversion process adopts a two-stage strategy: the first stage performs a coarse-grid search 502, traversing within the physically feasible region and locking the initial parameters 503; the second stage uses the initial parameters as a starting point, employs the LM Levenberg-Marquardt local optimization algorithm 504, combined with physical constraints / regularization 505, and guides the iteration direction through Jacobian calculation 508. When the convergence criterion 506 is met, the optimal parameters, including the zero-point offset of the time delay, are output 507, and uncertainty evaluation is completed simultaneously 509.

[0133] pass Figure 5 The process shown in this invention transforms the complex radar echo physical model ( Figure 4 This is embedded into a nonlinear inversion framework. By utilizing physical constraints and uncertainty assessment, not only is high-precision water level extraction achieved, but the technical challenges of traditional empirical heavy trackers being unable to handle complex inland water disturbances and lacking rigorous error estimation are also simultaneously solved.

[0134] S6. Distance and Elevation Calculation and Standard Correction

[0135] Based on the optimal water level time delay parameters, the slant distance and initial ellipsoidal height are calculated. The water surface ellipsoidal elevation is obtained by combining satellite attitude data. Corrections for wet / dry troposphere, ionosphere, solid tide, load tide, polar tide, and atmosphere are uniformly applied, and the water level is converted to a unified reference using a selected geoid model (such as EGM2008). This processing ensures the comparability of results from different orbits and time periods, laying the foundation for constructing long-term time series.

[0136] like Figure 6 As shown, the complete processing chain for converting raw time data from radar observations into absolute elevation data with hydrological significance is illustrated.

[0137] This process uses the round-trip time delay obtained from retracking as a basis, combines it with satellite orbit determination data for geometric calculation, and then progressively isolates the effects of atmospheric transmission delay and geophysical effects to ultimately achieve a unified vertical reference. Specifically, it includes the following key processing nodes:

[0138] (1) Basic geometric solutions (601, 602, 615, 603):

[0139] Round-trip delay 601: by Figure 5 The optimal time parameters are calculated using the parameter inversion steps shown. These are the raw observations of the radar pulses traveling back and forth between the satellite and the water surface.

[0140] Slant range calculation 602: Based on the speed of light constant and round-trip time delay, calculate the straight-line distance from the satellite center to the water surface reflection point.

[0141] Satellite Ellipsoid Elevation 615: The height of the satellite's center of mass relative to the reference ellipsoid, provided by precise orbit determination data, is the spatial reference point for calculating the Earth's surface elevation.

[0142] Water surface ellipsoid elevation 603 (uncorrected): Preliminary elevation value obtained by subtracting the slant distance from the satellite ellipsoid elevation. At this point, the data has not yet been adjusted for atmospheric delay and geophysical deformation, and only represents the relative position in a geometric sense.

[0143] (2) Geophysical corrections (604-610, 611, 612):

[0144] Various geophysical corrections (604-610): These include a series of error terms that must be subtracted to restore the true water level. Specifically, they include:

[0145] Atmospheric transport corrections: Moist troposphere 604, dry troposphere 605 and ionosphere 606 corrections are used to compensate for distance deviations caused by the reduced propagation speed of radar waves in the atmosphere.

[0146] Earth Deformation Correction: Corrections for Solid Tide 607, Load Tide 608, and Polar Tide 609 are used to eliminate the effects of vertical crustal movement caused by gravity.

[0147] Dynamic Atmospheric Correction 610: Used to eliminate the depressurization or elevation effect of atmospheric pressure changes on water surface height.

[0148] Correction summary node 611: Performs algebraic summation on all the above independent scalar correction items to generate the total correction.

[0149] Corrected ellipsoid elevation 612: By combining the uncorrected elevation with the corrected elevation, the pure water surface ellipsoid elevation after eliminating environmental interference is obtained.

[0150] (3) Reference conversion and final output (613, 614):

[0151] Geoid Model 613: Introduces a high-precision gravity field model to provide elevation anomalies between the reference ellipsoid and the geoid.

[0152] Unified benchmark water level 614: Final output product. By subtracting the geoid difference from the corrected ellipsoidal elevation, the elevation benchmark is transformed from a purely mathematically defined ellipsoid into a physically defined geoid (i.e., altitude), giving it practical value for cross-regional comparison and hydrological analysis.

[0153] The optimal round-trip time delay 601 obtained from the inversion is converted into the slant distance solution from the satellite to the water surface 602, initially obtaining the water surface ellipsoid elevation 603 (uncorrected). Subsequently, geophysical corrections such as moist tropospheric correction 604, dry tropospheric correction 605, ionospheric correction 606, as well as geophysical effect corrections 607, load tide correction 608, and polar tide correction 609 are applied item by item. Each item is superimposed at the correction summary node 611 to generate the corrected ellipsoid elevation 612, and combined with geoid models such as EGM2008 613, it is converted into a water level with a unified reference level 614.

[0154] pass Figure 6 The process shown in this invention enables a step-by-step transformation from the "time domain" to the "spatial geometric domain" and then to the "physical elevation domain." Through geophysical correction and benchmark unification, systematic errors introduced by path delay and crustal deformation are effectively eliminated, ensuring the consistency and accuracy of the final generated water level data under different spatiotemporal conditions.

[0155] S7. Multi-indicator quality control

[0156] Calculate the explanatory power of the simulated waveform for the observed waveform, and the whitening degree of the fitting residual; based on the explanatory power and the whitening degree, perform quality judgment and screening of the inversion results.

[0157] This invention employs a joint quality control system of "model explanation rate - residual whitening degree" for epoch-level result judgment and anomaly removal. The model explanation rate measures the ability of the echo numerical simulation to interpret observed energy; the residual whitening degree is used to verify whether the error is close to random noise rather than structural bias (such as significant autocorrelation, systematic underfitting / overfitting). A preferred judgment strategy combines threshold and weighted scoring: a model explanation rate ≥ 0.75 and a residual whitening degree score ≥ 0.6 indicate a pass; otherwise, it is marked as suspicious and reconstructed (e.g., adjusting the water body threshold or enabling a weak background term). Those still failing to meet the standards are removed or downweighted. This process effectively identifies problematic samples such as bridge obstruction, snow and ice cover, and extreme specular reflection, ensuring the consistency of time series and the stability of quality in large-scale production.

[0158] like Figure 7 The diagram illustrates a complete decision-making process for automatically screening and grading inversion results based on statistical characteristics. This process introduces two indicators: "model explanation rate" and "residual whitening degree," and evaluates the confidence level of the parameter inversion results through a combination of hard threshold gating and soft weighted scoring. Specifically, it includes the following key processing nodes:

[0159] (1) Input of indicators and initial screening of thresholds (701, 702, 708, 709):

[0160] Model Explanation Rate 701: One of the input metrics. It is used to quantify the goodness of fit between the numerical simulation waveform and the energy distribution of the real observed waveform. The closer the value is to a specific ratio (such as 1 or 100%), the higher the fidelity of the physical model to the surface scene.

[0161] Residual whitening degree 702: One of the input metrics. Used to test whether the fitted residual sequence conforms to a Gaussian white noise distribution. A high whitening degree indicates that the model has fully extracted the useful signal; a low whitening degree suggests the existence of systematic errors or interference that have not been captured by the model.

[0162] Explanation rate threshold gate 708 and whitening threshold gate 709: These are the "circuit breaker" thresholds set for the two indicators mentioned above. Any result falling below the preset physical reasonableness threshold will be marked as a potential anomaly to prevent severely distorted data from entering subsequent stages.

[0163] (2) Comprehensive scoring and decision-making (703, 704, 705):

[0164] Gated aggregation node 703: Logical connection point. It aggregates the indicator status after initial threshold screening with the original values, preparing for multi-dimensional fusion.

[0165] Weighted scorer 704: Constructs a comprehensive evaluation function based on the weights of each indicator's impact on water level accuracy. This step converts the interpretation rates of different physical dimensions and the degree of bleaching into a unified quality score.

[0166] Comprehensive Judgment 705: Final Decision Node (Rhombus). The calculated comprehensive score is compared with the quality pass line to determine whether the inversion result of the current epoch is usable.

[0167] (3) Output results and feedback (707, 706):

[0168] Quality Control Pass 707: Output status when the overall judgment is "Yes". This indicates that the water level data of this epoch is of reliable quality and can be archived as valid observations.

[0169] Reconstruction Triggered / Failed 706: Output status when the synthesis decision is "No". At this time, the data is marked as invalid. The system can choose to directly remove this point, or adjust the initial value based on the residual characteristics and trigger a new round of parameter inversion iteration (i.e., feedback to...). Figure 5 (A process) to attempt to correct the fitting bias.

[0170] pass Figure 7 The process shown in this invention establishes an automated quality control mechanism. By jointly evaluating the "overall similarity" (explanation rate) and "detail fidelity" (whitening degree) of waveform fitting, it effectively identifies and eliminates inferior data affected by complex terrain or strong interference, ensuring that the final generated water level product has high reliability and consistency.

[0171] S8. Product Output: Output a water level product that includes the absolute water level, inversion uncertainty, and quality indicator.

[0172] This invention outputs river water levels, epoch-level uncertainties, and quality indicators under a unified benchmark at a virtual station scale, and provides the energy proportion of "component waveforms—geographical sources" (main channel / floodplain / adjacent water bodies) to support physical interpretation. Water levels are standardized using unified geophysical corrections and elevation benchmarks to ensure comparability across orbits and time periods; uncertainties are estimated using a posteriori covariance and given in the form of one standard deviation; quality indicators are encoded using bitmasks to record whether each quality control indicator meets the standards and the type of anomaly. The product supports the fusion and offset correction of data from different orbits (such as based on cross-point constraints or reference water gauge calibration) to form continuous and stable river water level time series, meeting the needs of operational applications and model assimilation.

[0173] like Figure 8 As shown, Figure 8 This demonstrates the complete post-processing chain from discrete data from a single observation to the generation of high-quality, long-term hydrological products. The process aggregates scattered virtual station data, introduces multi-track fusion and external constraint mechanisms to eliminate system bias, and after time-series optimization, finally generates standardized water level time-series products. Specifically, it includes the following key processing nodes:

[0174] (1) Data aggregation and access (801, 803, 804, 800):

[0175] Virtual station water level 801, epoch-level uncertainty 803, and quality mark 804: these are the underlying input units for product generation. Among them, the quality mark 804 typically uses a bitmask to record the status information of the data at each stage of inversion and quality control; uncertainty 803 provides error statistics for subsequent weighted fusion.

[0176] Product output summary node: Logical aggregation point. Responsible for summarizing the outputs from previous steps (such as...). Figure 7 (Quality control) Select qualified single-epoch observation data and collect them according to spatial location under the corresponding virtual station name.

[0177] (2) Multi-source fusion and bias correction (805, 807, 808, 806):

[0178] Multi-track fusion 805: For situations where the same body of water may be covered by multiple satellite orbits (within the same orbit, in different orbits, or by different satellite missions), observation sequences from different orbits are spliced ​​together in time and space.

[0179] Reference water gauge / control point 807 and intersection constraint 808: external and internal constraints. Using measured water gauge data or water level differences at satellite orbit intersections, calculate systematic deviations between different orbits or missions.

[0180] Offset / Scale Correction 806: Based on the correction parameters calculated under the above constraints, systematic errors are eliminated from the original water level sequence to ensure the compatibility of multi-source data under the same benchmark.

[0181] (3) Timing optimization and product generation (809, 812, 802):

[0182] Time Series Consistency Screening / Outlier Removal 809: A secondary quality control measure performed on the time series dimension. It identifies and removes outliers that statistically deviate significantly from historical trends, further improving the smoothness and reliability of the series.

[0183] Time resampling / interpolation 812: Optional processing step. Resamples non-uniformly sampled observation data to regular time intervals to facilitate hydrological analysis or model assimilation.

[0184] Water Level Time Series 802: The final unified benchmark core product, which includes historical water level change data that has been corrected, cleaned and standardized.

[0185] (4) Management and Release (811, 810):

[0186] Metadata and Version Management 811: Attach descriptive information to data products (such as processing algorithm version, coordinate system definition, data coverage) and establish a version tracking mechanism to support traceability.

[0187] Product Archiving and Release 810: Process Endpoint. The packaged standard data product is stored in the database and released to users through the service interface.

[0188] The calculation of indicators such as model explanation rate (701) and residual whitening degree (702) is quantified by the weighted scorer (704) and then enters the comprehensive judgment (705). If the judgment result is quality control passed (707), the data flows into the product output summary node (800) to generate the final water level time series (802), epoch-level uncertainty (803), and quality indicator (804). If the judgment result is reconstruction triggered / failed (706), the data is fed back to the scenario construction module (300) to adjust parameters and re-invert. Finally, the qualified product undergoes multi-track fusion (805) and time series consistency screening (809) before being archived and released (810).

[0189] pass Figure 8The process illustrated demonstrates how this invention bridges the gap between "instantaneous observation" and "long-term time-series monitoring." Through multi-track fusion and a rigorous deviation correction mechanism, it effectively solves the problem of inconsistent benchmarks for multi-source satellite data, generating high-temporal-resolution and high-precision time-series products of inland water levels.

[0190] To verify the feasibility and technical effectiveness of this invention, four representative embodiments and comparative results are provided. The embodiments are consistent with the aforementioned technical solutions, covering the extreme and intermediate values ​​of key hyperparameters, and performance characterization indicators are given. Unless otherwise specified, all processing procedures follow... Figure 1 The architecture shown is executed.

[0191] 1. General Test Environment and Setup

[0192] (1) Data and Task: 20 Hz waveform data in Sentinel-3 / SRAL SAR mode, according to... Figure 1 Each functional module executes a complete processing chain.

[0193] (2) Software and hardware: 64-bit CPU (≥8 cores), memory ≥16 GB; parallel acceleration is enabled for grid integration and Jacobi calculation; the optimizer is Levenberg-Marquardt (LM).

[0194] (3) Geometric and sampling parameters: 256 distance gates for unified time delay grid; the footprint grid resolution is adaptively selected at three levels of 5 meters, 15 meters and 30 meters according to the river width, corresponding to narrow, medium and wide river scenarios respectively.

[0195] (4) Water body identification parameters: The frequency thresholds of water bodies appearing in GSWO were tested at 0.10, 0.20, and 0.35 (covering the end value and median value); only water bodies were modeled, and non-water bodies were not included in the integration; if necessary, the weak background term A_bg∈{0%,3%} was enabled to absorb residual non-water energy.

[0196] (5) Inversion and constraints: The coarse search range of water level is set to ±2 meters, ±3 meters, ±5 meters, with step sizes of 0.05 meters, 0.10 meters, and 0.20 meters; the maximum iteration of LM is 50, and the convergence threshold is 1e-6; physical constraints include water level-DEM consistency and roughness feasible range.

[0197] (6) Quality control thresholds: Model explanatory power ≥ 0.75; residual whitening ≥ 0.6; a quality control strategy combining threshold judgment and weighted scoring is adopted (see attached table). Figure 7 (The process).

[0198] Example 1: Benchmark River Section (Simple River, Simple Scenario)

[0199] This embodiment selects a typical straight river section (geographically approximately 44.42°N–44.47°N, 0.16°E–0.22°E) near the Marmande hydrological station on the Garonne River in France as the benchmark control area. This river section is approximately 150–250 m wide, and there is no significant interference from nearby still water bodies or wide floodplains within the effective footprint of the radar altimeter. The electromagnetic scattering environment is stable, making it suitable for verifying the basic performance of the method of this invention under low-interference conditions. According to the method described in this invention, the experimental parameters are set as follows: the footprint is discretized using a 5 m grid resolution; a water frequency threshold of O=0.10 is set using the GSWO global surface water dataset to extract the mask; in the parameter inversion stage, the coarse search range is set to ±3 m, and the search step size is 0.10 m.

[0200] During the evaluation period in 2024, the consistency between the simulated waveforms generated by this invention and the satellite-measured waveforms will first be quantitatively evaluated. For example... Figure 9 The figure shows a comparison between the numerically simulated waveform and the satellite-measured waveform in a typical river section. The figure uses normalized echo power as the vertical axis and the range gate (or sampling time) as the horizontal axis, visually demonstrating the degree of fit between the model's forward modeling results and real radar observation data. The figure typically contains two main curves: the solid blue line represents the measured waveform received by a satellite altimeter (such as Sentinel-3 SRAL), and the dashed red line represents the optimal simulated waveform generated based on the physical model of this invention. Statistical results show that the average correlation coefficient (r) of the waveform fitting reaches 0.89, the Kling-Gupta efficiency coefficient (KGE) is 0.71, and the normalized root mean square error (NRMSE) is 0.04. These indicators collectively demonstrate that the numerically simulated waveform constructed in this invention exhibits good physical consistency with the measured data in terms of morphology, energy amplitude, and variation trends, verifying the applicability of the forward modeling model in simple river scenarios. To evaluate the accuracy of water level inversion, ground-measured water gauge data were used as the true reference. The method of this invention was compared with two baseline retracking algorithms, OCOG (centroid shift method) and SAMOSA (physical model method), under the same data caliber. Evaluation metrics included root mean square error (RMSE), normalized root mean square error (NRMSE), mean absolute error (MAE), correlation coefficient (r), and KGE coefficient. Figure 10 As shown, Figure 10 This is a time series comparison of water level inversion results between the method of this invention and existing mainstream retracking algorithms. The figure uses time (year / month) as the horizontal axis and absolute water level elevation based on the reference ellipsoid as the vertical axis, showing the long-term monitoring performance of different algorithms at the same hydrological station.

[0201] Table 1 presents representative comparison results. Overall, the inversion results of each algorithm are basically consistent with the measured water level change trends. Regarding accuracy, the method of this invention exhibits excellent stability: the correlation coefficient r of this invention is 0.93, slightly better than the OCOG algorithm (0.92) and the SAMOSA algorithm (0.92); the normalized root mean square error (NRMSE) of this invention is 0.26, slightly lower than the OCOG algorithm (0.27) and the SAMOSA algorithm (0.27).

[0202] Table 1 Comparison of water level inversion accuracy of different algorithms in the benchmark river section

[0203]

[0204] Conclusion: Experimental results show that in simple scenarios involving only rivers, the method of this invention performs comparably to the classic OCOG and SAMOSA algorithms, effectively achieving accurate water level inversion. This demonstrates the correctness of the underlying physical model of this invention, and its ability to degenerate into a standard model without complex disturbances, exhibiting good benchmark compatibility and inversion reliability. This provides a reliable algorithmic baseline for subsequent processing of complex scenarios.

[0205] Example 2: Simulation Experiment of Waveform Response Mechanism in Typical Multi-Water Body Scenarios

[0206] This embodiment constructs four typical radar altimetry footprint scenarios (corresponding to the attached diagrams). Figure 11 By using the controlled variable method, the specific effects of physical factors such as water volume, spatial distance, elevation difference, and effective area on the echo waveform were simulated and analyzed to verify the numerical simulation method of the present invention's ability to finely characterize complex scenarios and its physical interpretability.

[0207] Figure 11The figures show the numerical simulation waveforms for four typical multi-water-body interference scenarios. These figures correspond to the mechanism verification experiments conducted using the controlled variable method in the embodiments, and include four sub-figures, each illustrating the physical model's response characteristics to different surface geometric elements: Sub-figure (a) Baseline Single Water Body Scenario: Shows the standard specular reflection waveform when only a single water body at the nadir point is included, serving as a reference benchmark in the interference-free state. Sub-figure (b) Distance Factor Response Scenario: Shows the waveform after introducing an interfering water body at the same height as the nadir point but with a horizontal offset (e.g., 500 m). A significant second sub-peak is visible, with a delayed peak position and amplitude attenuation, reflecting the combined effect of off-axis distance on echo delay and energy. Sub-figure (c) Elevation Difference Response Scenario: Shows the waveform when there is a vertical elevation difference (e.g., 5.0 m) between the main water body and the interfering water body. The two sub-peaks show significant misalignment on the distance gate axis, with a step at the waveform leading edge, demonstrating the decisive role of elevation difference in echo peak separation. Sub-figure (d) Area difference response scenario: This shows the waveform when the effective area of ​​the interfering water body is doubled (e.g., 2 times). It can be seen that the relative amplitude ratio of the two sub-peaks changes, reflecting the linear contribution of the effective scattering area to the echo energy.

[0208] In the baseline single water body scenario ( Figure 11 In a), it is assumed that there is only a single water body at the nadir point within the radar footprint. The normalized power-range gate curve obtained by simulation exhibits the quasi-specular reflection characteristics of a standard inland water body. Its steep leading edge and sharp main peak provide a reference benchmark for subsequent multi-water body interference analysis.

[0209] To address the complex situation of two coexisting water bodies, the experiment investigated the specific effects of distance, elevation, and area factors on waveform morphology. For the two-water-body-distance-factor scenario (… Figure 11 (b) After introducing a body of water at the same height as the nadir point but about 500 meters away horizontally, the simulated waveform exhibits a distinct double-peak structure; the second sub-peak is characterized by a lagging position and a lower amplitude due to the increased path delay caused by off-axis observation geometry and the attenuation of antenna gain.

[0210] Based on this, the scenario of two water bodies—elevation difference ( Figure 11 c) Keeping the horizontal geometric parameters unchanged, but setting a 5.0-meter elevation difference between the main water body and the interfering water body, the results show that the two sub-peaks are significantly misaligned and separated on the distance gate axis, accompanied by a step characteristic of the waveform leading edge, which intuitively reflects the decisive role of the vertical elevation difference on the echo delay.

[0211] For scenarios with differences in effective incident area ( Figure 11 d) Under the conditions of equal height and equal distance, the area of ​​the interfering water body is set to twice that of the target water body. The simulation results show that the relative amplitude ratio between the two changes significantly, and the energy of the sub-peak corresponding to the interfering water body is significantly enhanced, reflecting the linear contribution of the effective reflection area to the backscattered energy.

[0212] In summary, simulation experiments have demonstrated that the component echo model constructed in this invention has clear physical interpretability and can effectively distinguish and quantify the specific impact of physical factors such as distance, elevation, and area on signal aliasing in multi-water scenarios. This provides a solid physical model foundation for solving waveform re-tracking problems under conditions of interference from adjacent lakes or complex tidal flats.

[0213] A river level inversion system that takes into account the heterogeneity of altimeter footprints, such as Figure 1 As shown, it includes:

[0214] The data access and preprocessing module 100 is used to acquire and preprocess the observation waveform data, satellite orbit parameters, attitude parameters, instrument parameters, geophysical correction data and space prior data of the synthetic aperture radar altimeter;

[0215] The geometric modeling module 200 is used to discretize the radar ground illumination area into a grid and calculate the geometric parameters of each grid cell.

[0216] Scene construction module 300 is used to identify and classify water body pixels within the footprint grid, mainly river water bodies, floodplain / tidal flat water bodies, and adjacent isolated water bodies;

[0217] The numerical echo simulation module 400 is used to perform electromagnetic scattering numerical integration on the footprint grid based on differentiated surface roughness parameters and geometric parameters to generate simulated echo waveforms.

[0218] The parameter inversion module 500 is used to solve for the optimal water level delay parameters that minimize the difference between the simulated waveform and the observed waveform through coarse grid search and local optimization.

[0219] The solution and correction module 600 is used to calculate the absolute water level based on the optimal parameters and perform geophysical corrections.

[0220] The quality control module 700 is used to judge and screen the inversion results based on the interpretation rate of the simulated waveform and the whitening degree of the fitting residual.

[0221] Product output module 800 is used to output water level products that include absolute water level, inversion uncertainty and quality indicators.

[0222] An electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the steps of the above-described method for river level inversion that takes into account the heterogeneity of altimeter footprints.

[0223] In the above embodiments, the descriptions of each embodiment have different focuses. For parts that are not described in detail or recorded in a certain embodiment, please refer to the relevant descriptions of other embodiments.

[0224] The above-described embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.

Claims

1. A method for river water level inversion that takes into account the heterogeneity of altimeter footprints, characterized in that, Includes the following steps: S1. Data Acquisition and Preprocessing: Acquire observation waveform data from synthetic aperture radar altimeter, satellite orbit parameters, attitude parameters, instrument parameters, geophysical correction data, and spatial prior data including river network, water body mask, and digital elevation model; The observed waveforms are subjected to noise removal and normalization, and then mapped to a uniform time delay grid. S2. Footprint Gridding and Observation Geometric Modeling: The radar ground illumination area corresponding to a single observation epoch is discretized into a footprint grid; based on the satellite orbital parameters, attitude parameters, and the digital elevation model, the geometric parameters of each grid cell relative to the satellite are calculated, including slant range, incident angle, and azimuth angle; S3. Water body scene construction: Within the footprint network determined in step S2, water body pixels are identified based on the water body mask; according to the connectivity of the water body pixels with the centerline of the river network and the spatial position of the water body pixels on the digital elevation model, the water body pixels are classified. The classification results of the water body pixels are: main river channel water body, floodplain / tidal flat water body connected to the main river channel, and adjacent isolated water body not connected to the main river channel. S4. Numerical echo simulation with electromagnetic scattering constraints: Different surface roughness parameters are set for the main channel water body and adjacent isolated water bodies respectively; Based on the geometric parameters calculated in step S2 and the classification results in step S3, electromagnetic scattering numerical integration is performed on the footprint grid to generate simulated echo waveforms; wherein, water level changes are simulated by changing the elevation of the grid cells corresponding to the main channel water body, and the elevation changes are converted into the overall translation of the simulated waveform on the time delay axis through geometric relationships; S5. Waveform re-tracking and parameter inversion: Construct an objective function that minimizes the difference between the simulated waveform and the observed waveform; perform a coarse grid search in the interval near the initial time delay value determined based on the digital elevation model to obtain the initial value of the inversion parameters; perform local optimization solution starting from the initial value to obtain the optimal water level time delay parameters; S6. Distance and Elevation Calculation and Standard Correction: Based on the optimal water level time delay parameters, calculate the slant distance and initial ellipsoidal height, and apply wet / dry tropospheric correction, ionospheric correction, and solid tide correction one by one, and convert to the absolute water level under the unified geoid reference; S7. Multi-index quality control: Calculate the explanatory power of the simulated waveform to the observed waveform, and the whitening degree of the fitting residual; perform quality judgment and screening on the inversion results based on the explanatory power and the whitening degree; S8. Product Output: Output a water level product that includes the absolute water level, inversion uncertainty, and quality indicator.

2. The method according to claim 1, characterized in that, In step S3, the classification of water body pixels is specifically as follows: water body pixels that are connected to the center line of the river network and located within the main channel space are identified as main river channel water bodies; water body pixels that are connected to the main river channel water bodies but located outside the main channel space are identified as floodplain or tidal flat water bodies; and water body pixels that exist within the footprint grid but have no topological connection with the main river channel water bodies are identified as adjacent isolated water bodies.

3. The method according to claim 1, characterized in that, In step S4, the surface roughness parameter is a mean square slope parameter based on a geometric optical approximation model, used to characterize the statistical properties of the microscopic undulations of the water surface; wherein, "the mean square slope parameter value set for adjacent isolated water bodies" is less than "the mean square slope parameter value set for the main river channel water body".

4. The method according to claim 1, characterized in that, The electromagnetic scattering numerical integration in step S4 specifically includes: for each water body grid cell, calculating its backscattering coefficient according to the corresponding surface roughness parameter based on the category, and calculating the range attenuation factor and antenna pattern gain weight in the time delay-Doppler domain in combination with geometric parameters, and generating the component waveform of this type of water body through numerical integration; and performing incoherent linear superposition of the component waveforms of the main channel water body, the adjacent isolated water body, and a weak background noise term to generate the simulated echo waveform.

5. The method according to claim 1, characterized in that, In step S5, the range of the coarse grid search is determined based on the prior water surface elevation provided by the digital elevation model, combined with the time delay reference value obtained by satellite orbital altitude conversion, and fluctuates up and down by a preset range of physical water level changes, which is 2 meters to 10 meters.

6. The method according to claim 1, characterized in that, In step S5, the local optimization solution adopts the Levenberg-Marquardt algorithm, and during the optimization process, a value range constraint based on physical experience is applied to the surface roughness parameter to prevent the parameter optimization from falling into non-physical solutions.

7. The method according to claim 1, characterized in that, Step S5 further includes: calculating and outputting the inversion uncertainty corresponding to the water level delay parameter based on the posterior covariance information obtained from the optimization process.

8. The method according to claim 1, characterized in that, In step S7, the quality determination is as follows: when the explanatory rate is not lower than the first threshold and the whitening degree is not lower than the second threshold, the epoch inversion result is determined to be of qualified quality; otherwise, it is discarded or marked. The first threshold is 0.75 and the second threshold is 0.

6.

9. The method according to claim 1, characterized in that, Step S8 also includes: fusing multiple qualified water level products obtained by inverting the same virtual station from different satellite orbits, and using orbital intersection consistency constraints or ground reference water level data to perform system deviation correction, thereby generating a time-series continuous river water level product.

10. A river level inversion system that takes into account the heterogeneity of altimeter footprints for performing the method of claim 1, characterized in that, include: The data access and preprocessing module is used to acquire and preprocess the observation waveform data, satellite orbit parameters, attitude parameters, instrument parameters, geophysical correction data and space prior data of the synthetic aperture radar altimeter; The geometric modeling module is used to discretize the radar ground illumination area into a grid and calculate the geometric parameters of each grid cell. The scene construction module is used to identify and classify water body pixels within the footprint grid, mainly river water bodies, floodplain / tidal flat water bodies, and adjacent isolated water bodies. The numerical echo simulation module is used to perform numerical integration of electromagnetic scattering on the footprint grid based on differentiated surface roughness parameters and geometric parameters, and generate simulated echo waveforms. The parameter inversion module is used to solve for the optimal water level delay parameters that minimize the difference between the simulated waveform and the observed waveform through coarse grid search and local optimization. The solution and correction module is used to calculate the absolute water level based on the optimal parameters and perform geophysical corrections. The quality control module is used to judge and screen the inversion results based on the interpretation rate of the simulated waveform and the whitening degree of the fitting residuals. The product output module is used to output water level products that include absolute water level, inversion uncertainty, and quality indicators.

Images