Underwater terrain inversion method based on height measurement satellite

Abstract

The invention provides an underwater terrain inversion method based on a height measurement satellite, belongs to the technical field of remote sensing geoscience application, and constructs an automatic processing flow integrating space screening, elevation correction, signal denoising, double-strategy surface and bottom recognition and physical correction. The method is suitable for underwater topography generation in complex water body environments such as lakes and reservoirs, and the efficiency of underwater topography data acquisition is greatly improved. A strategy of combining adaptive kernel density estimation and clustering analysis is adopted to remarkably improve the water surface and water bottom signal recognition capability, and a multi-temporal collaborative analysis and refractive index physical correction model is introduced to realize high-precision underwater terrain extraction of water bodies such as inland reservoirs. The method has obvious advantages in the aspects of automation degree, application range and result stability, and has important engineering application and popularization value.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of remote sensing geoscience application, and in particular to an underwater topography inversion method based on an altimetry satellite. BACKGROUND

[0002] Underwater topography information is an important basic data for water resource management, flood control and disaster prevention, ecological protection and engineering planning. Existing underwater topography measurement mainly relies on shipborne multi-beam sonar systems or manual measurement methods, which have high accuracy but generally have problems such as large workload, high cost, limited coverage, etc., and are difficult to meet the rapid monitoring needs of large-scale water areas. In recent years, satellite remote sensing technology has provided a new solution for bathymetry: optical remote sensing can estimate shallow water depth by water color inversion, and laser depth sounding satellites can directly detect water surface and water bottom reflection signals to realize underwater topography inversion. However, in complex water areas such as inland lakes and reservoirs, existing methods still face significant challenges: (1) optical images are easily affected by water turbidity and suspended matter, resulting in a significant decrease in inversion accuracy; (2) laser depth sounding signals are easily disturbed by water surface reflection and high noise background, and water surface and water bottom reflections are difficult to separate stably, resulting in large calculation errors; (3) orbital photon data processing generally relies on manual interpretation, has low automation, lacks geophysical correction, and is difficult to adapt to complex terrain and multi-temporal data analysis needs. SUMMARY

[0003] The purpose of the present application is to provide an underwater topography automatic inversion method based on the fusion of ICESat-2 ATL03 photon signals and high-resolution optical images, which solves the technical problems of large photon data noise interference, difficult surface-bottom identification and high manual involvement in the prior art.

[0004] An automatic processing flow is constructed, which integrates spatial screening, height correction, signal denoising, double-strategy surface-bottom identification and physical correction. It is suitable for underwater topography generation in complex water environments such as lakes and reservoirs, greatly improving the efficiency of underwater topography data acquisition.

[0005] The strategy of combining adaptive kernel density estimation with clustering analysis significantly improves the ability to identify water surface and water bottom signals, and introduces multi-temporal collaborative analysis and refractive index physical correction model to realize high-precision underwater topography extraction of inland reservoirs and other water bodies. The present application has obvious advantages in automation, application range and result stability, and has important engineering application and promotion value.

[0006] In order to achieve the above-mentioned purpose, the technical scheme adopted by the present application is as follows:

[0007] An underwater topography inversion method based on an altimetry satellite, the method comprising the following steps:

[0008] Step 1: Extract the water body vector boundary from the optical image data within the corresponding time range;

[0009] Step 2: Select ICESat-2 ATL03 data covering the water body, read the photon longitude, latitude, ellipsoidal height, signal confidence parameters of the specified beam from the data, and perform spatial filtering combined with the vector boundary;

[0010] Step 3: Obtain the geodetic datum anomaly value in the geophysical correction parameters corresponding to the photon, correct the photon height through the geodetic datum anomaly value, and match it to each photon point through the interpolation method to realize the conversion from ellipsoidal height to orthometric height;

[0011] Step 4: Photon signal processing stage, through quantile filtering and denoising, retain weak signals with water bottom greater than a certain value and suppress upper noise point interference;

[0012] Step 5: Separate the water surface and underwater photons, fuse the density-based spatial clustering algorithm and the double-strategy bottom identification mechanism based on kernel density estimation in the sliding window, use DBSCAN clustering to separate the water surface cluster and the water bottom cluster, and when the DBSCAN clustering result is unstable or fails, automatically switch to the KDE strategy to extract the double peak value, ensuring the robustness and stability of the surface-bottom identification;

[0013] Step 6: Establish a refractive index physical correction model based on the difference in light propagation speed between water and air, and perform refractive correction on the original water depth data to obtain the true water bottom height;

[0014] Step 7: Perform median filtering and smoothing processing on the water bottom height sequence of the entire orbit, and output several-dimensional parameter data sets including longitude, latitude, water surface height, original water bottom height, and corrected water bottom height;

[0015] Step 8: Integrate several orbit data, fuse and exclude outliers for all transit orbit data within a fixed time period, and generate continuous underwater terrain based on interpolation;

[0016] Step 9: Use several beam sonar measured underwater terrain data to verify the obtained underwater terrain, compare the corrected water bottom height with the measured data, and evaluate the inversion accuracy.

[0017] Further, in step 1, the time of optical image data and the time of altimetry satellite data are controlled within 1 month, and the optical image data includes Landsat and Sentinel-2 optical satellite data.

[0018] Further, in step 1, the normalized water body index optical image is used to extract the water body vector boundary, and the wavebands used in different satellite data sources have some differences. The normalized water body index NDWI formula is as follows:

[0019]

[0020] wherein represents the green wave band, corresponding to the B2 band of Landsat 5 and Landsat 7, corresponding to the B3 band of Landsat 8, Landsat 9 and Sentinel-2, represents the near-infrared wave band, corresponding to the B4 band of Landsat 5 and Landsat 7, corresponding to the B5 band of Landsat 8 and Landsat 9, corresponding to the B8 band of Sentinel-2.

[0021] Further, in step 2, the ICESat-2 ATL03 data is in h5 format, the strong beam photon information is selected, the photon points in the water body range are screened in combination with the water body vector boundary, and the output is in Shapefile format, containing the parameters of longitude, latitude, ellipsoidal height, and signal confidence.

[0022] Further, in step 4, in the photon signal processing stage, since the water surface photons are greater than the set value and easy to identify, a random denoising strategy based on height distribution is adopted with a sliding window of 75 photons, the upper noise is reduced according to the height distribution, and the weak signal at the bottom is highlighted.

[0023] Further, in step 6, the specific process of establishing the light propagation refractive index correction model is to correct the preliminary water depth result according to the difference between the air and water propagation speed, to obtain the true water bottom elevation, and the refractive correction model is corrected based on the ratio of the air refractive index to the water refractive index , as follows:

[0024]

[0025] wherein is the corrected water bottom elevation, wherein is the water surface elevation, is the original water bottom elevation.

[0026] Further, in step 7, the water bottom elevation sequence of the whole orbit is subjected to median filtering and smoothing processing to eliminate local abnormal points, forming a continuous and stable underwater topographic profile, and the data set is used for underwater topographic modeling and mapping.

[0027] Further, in step 8, all transit orbit data from 2019 to 2025 are integrated, and abnormal underwater elevations are removed, and then a digital elevation model is generated based on an interpolation method.

[0028] Further, in step 9, the inversion result is verified in precision by using the underwater topographic data measured by the multi-beam sonar, and in the verification process, the correction of the geodetic datum anomaly between data needs to be paid attention to, the inversion precision is evaluated by comparing the corrected water bottom elevation with the sonar measured result, so that the reliability of the method is ensured.

[0029] The present application has the following beneficial effects due to the adoption of the above technical solutions:

[0030] (1) The present application has significant advantages in automation degree, precision stability, application range and multi-time phase adaptation capability, and is particularly suitable for large-scale and rapid mapping tasks in complex inland waters.

[0031] (2) The present application has important engineering popularization value and intellectual property protection prospect, and can be widely applied in the fields of water resource management, channel design, reservoir scheduling, ecological monitoring, disaster prevention and reduction, and digital twin river basin construction.

[0032] (3) Economic benefits: compared with the traditional shipborne measurement method, the present application can reduce the measurement cost per square kilometer by more than 70%, and save a large amount of labor and equipment cost in large area operation.

[0033] (4) Social benefits: the present application can improve the efficiency of water resource management, reduce the interference of field operation on ecological environment, and improve the water disaster warning and emergency response capability. BRIEF DESCRIPTION OF DRAWINGS

[0034] Figure 1 is the overall flow chart of the underwater topographic inversion method of the present application;

[0035] Figure 2 is the schematic diagram of the optical image extraction water boundary and photon screening of the present application;

[0036] Figure 3 is the schematic diagram of photon denoising of the present application;

[0037] Figure 4 is the schematic diagram of the content contained in the output data of the present application;

[0038] Figure 5 is the schematic diagram of the water bottom height sequence smoothing and DEM generation of the present application;

[0039] Figure 6 is the comparison and verification diagram of the inversion result and multi-beam measured data of the present application. DETAILED DESCRIPTION

[0040] For the purposes of the present application, the technical solutions and advantages will be more clearly apparent from the following detailed description of preferred embodiments, given by way of example and with reference to the accompanying drawings. However, it should be noted that the numerous details set forth in the specification are merely provided to give a thorough understanding of one or more aspects of the present application and are not intended to limit the present application to these aspects.

[0041] The example data used in the present application includes Sentinel-2 MultiSpectral Instrument Level-2A (SR) image data provided by the European Space Agency with a spatial resolution of 10x10 m, ATL03 products of ICESat-2 altimetry satellite provided by the U.S. Snow and Ice Data Center, and underwater topographic survey data from T50-P multi-beam bathymeter. The example Sentinel-2 image data, ICESat-2 altimetry data, and multi-beam underwater topographic survey data cover the Pingheshu Reservoir in Mugeng Town, Guiping City, Guigang City, Guangxi Zhuang Autonomous Region, and the time span is from January 1, 2018 to March 31, 2025. The underwater topographic survey data was provided by the Guangxi Zhuang Autonomous Region Water Resources and Electric Power Survey and Design Research Institute Co., Ltd., and was obtained in June 2025.

[0042] As shown in Figure 1 , the specific steps of the method include the following:

[0043] Step 1: Select Sentinel-2 image data as close as possible to the same month as ICESat-2 time, if there is no image of the same month, use the image data of the adjacent months instead, a total of 14 images are collected.

[0044] Step 2: Calculate the NDWI of the Sentinel-2 image, and use the OTSU best threshold method to segment the NDWI to obtain rough water body and non-water body raster information.

[0045] Step 3: Use ArcGIS to perform vectorization operation on the water body information, convert tiff format to shapefile format, remove vector areas smaller than three pixels, and correct the vector boundary in combination with visual interpretation to obtain complete and accurate water body vector boundary Figure 2 ).

[0046] Step 4: Reading and filtering ICESat-2 ATL03 data based on Python. Based on the h5py library, the photon latitude, longitude, ellipsoidal height, signal confidence, and geoid anomaly values in the ATL03 HDF5 file are read. To ensure data quality, this implementation prioritizes strong-beam photons and removes photons with a signal confidence of less than 2. After filtering and leveling, the photon data is output in CSV and Shapefile formats for subsequent analysis.

[0047] Step 5: Sort the filtered ATL03 photon point data by satellite orbit direction (ascending latitude), and use a sliding window analysis strategy with a sliding window of 75 photons (i.e., 75 photons before and after the current point, with a total window of 151 photons). Water surface and water bottom identification is performed separately within each window to ensure local stability of the results.

[0048] Step 6: Perform quantile denoising on the photon heights within the window. This method retains weak signal photons at the bottom and randomly removes non-bottom points at the top, reducing noise interference. Parameters are set to bottom percentile = 5% and reduction ratio = 0.5. As shown in Figure 3 , determine noise and valid photons.

[0049] Step 7: Use a two-strategy approach to identify the bottom, with DBSCAN clustering method as the first choice: use photon height as a one-dimensional feature for density clustering, automatically find high-density clusters corresponding to water surface and water bottom clusters, and use cluster median as the representative value of water surface and water bottom height. In the DBSCAN clustering method: input is the denoised photon height within the window; automatically search for the optimal clustering parameters eps ∈ [0.5,10], min_samples ∈ [2,5]; output two clusters corresponding to water surface and water bottom; if clustering fails or the number of clusters is insufficient, automatically switch to the KDE method.

[0050] Step 8: KDE method is used to extract the double peak value of height distribution, which is only used when DBSCAN cannot identify valid clusters to ensure the robustness of water surface and water bottom height identification. Gaussian kernel density estimation is used, with bandwidth automatically calculated as 30% of the standard deviation of photon height and limited to [0.1,1.0]; evaluation resolution = 0.01 m, max_points = 10000.

[0051] Step 9: Refraction correction, based on the physical law of light propagation between air (refractive index n1≈1.0) and water (refractive index n2≈1.333), the original water depth is corrected using the refractive index correction model proposed by Parrish. The corrected water bottom height is more in line with the actual situation, effectively reducing the system error caused by light refraction.

[0052] Step 10: If the calculated water depth is outside the reasonable range (for example, less than 0 m or greater than 60 m), it is determined to be abnormal and is rejected. The identified water depth sequence and water bottom height sequence are median filtered and smoothed (kernel_size = 5), and the final output results include: water surface height, water bottom height, corrected water bottom height, original water depth, and smoothed corrected water depth, etc.

[0053] The output data structure is as shown in Figure 4 Each photon point output field includes:

[0054] latitude: latitude;

[0055] longitude: longitude;

[0056] surface: water surface height;

[0057] bottom: preliminary water bottom height;

[0058] cor_bottom: refractive index corrected water bottom height;

[0059] depth_raw: original water depth;

[0060] cor_depth: smoothed and corrected water depth.

[0061] Step 11: Loop processing of photon CSV files for the same track or multiple tracks, automatically complete sorting, sliding window analysis, denoising, DBSCAN / KDE identification, refractive correction and smoothing output. Figure 5 Display the smoothed results. During batch processing, the input directory and output directory can be set, supporting parallel or serial processing of multiple files. The final interpolation obtains the underwater DEM result, as shown in Figure 5 The local terrain is shown.

[0062] Key parameter settings and optimization

[0063] Sliding window size window_size = 75 photons;

[0064] Quantile denoising bottom percentile = 5%, upper point deletion ratio reduction_ratio = 0.5;

[0065] DBSCAN eps ∈ [0.5,10], min_samples ∈ [2,5];

[0066] KDE bandwidth adaptive, limited in [0.1,1.0];

[0067] KDE evaluation resolution resolution = 0.01 m, maximum number of points max_points = 10000;

[0068] Refractive index n1 / n2 = 1.0 / 1.333;

[0069] Median filter kernel_size = 5;

[0070] Depth anomaly range max_depth = 60 m.

[0071] After the above process, the output underwater terrain data is smooth, continuous and the noise is significantly reduced. The underwater terrain accuracy is related to the water depth. When the estimated water depth is within 30 meters, 500 sample points are randomly selected, and the average deviation between the estimated water bottom elevation and the measured data of the multi-beam sonar is 5.27 m, as shown in the following table. Figure 6 This method realizes a great degree of automatic processing, reduces manual intervention, improves the efficiency and accuracy of underwater terrain inversion, and can be widely applied to the generation of high-precision underwater terrain of lakes, reservoirs and inland water bodies.

[0072] The parameters of each step can be adjusted according to different research areas and photon data characteristics, including sliding window size, MI denoising ratio, DBSCAN parameters and KDE bandwidth, etc.

[0073] The adaptive kernel density estimation and clustering analysis strategy significantly improves the signal recognition ability of water surface and water bottom, and introduces multi-temporal collaborative analysis and refractive index physical correction model to realize high-precision underwater terrain extraction of inland reservoirs and other water bodies. Compared with the prior art, the present application has obvious advantages in automation degree, application range and result stability, and has important engineering application and popularization value.

[0074] The remaining matters of the present application are known technologies.

[0075] The above description is only the preferred embodiments of the present application, and it should be pointed out that for ordinary skilled in the art, without departing from the principles of the present application, a number of improvements and refinements can be made, and these improvements and refinements should be considered as the protection scope of the present application.

Claims

1. A method for underwater topography inversion based on altimetry satellites, characterized in that, The method includes the following steps: Step 1: Select optical images within the corresponding time range to extract the water body vector boundary; Step 2: Select ICESat-2 ATL03 data covering the water body, and read the parameters of photon latitude and longitude, ellipsoidal height, and signal confidence of the specified beam from the data, and perform spatial filtering in combination with vector boundaries; Step 3: Obtain the geoid anomaly value in the geophysical correction parameters corresponding to the photon, correct the photon elevation using the geoid anomaly value, and match it to each photon point using an interpolation method to achieve the conversion of ellipsoidal height to orthographic height; Step 4: In the photon signal processing stage, weak signals above the set value at the bottom of the water are retained and interference from noise points at the top are suppressed through quantile filtering and noise reduction. Step 5: Perform surface and underwater photon separation. Within a sliding window, fuse density-based spatial clustering algorithm and kernel density estimation dual-strategy bottom identification mechanism. Use DBSCAN clustering to separate surface clusters and bottom clusters. When DBSCAN clustering results are unstable or fail, automatically switch to KDE strategy to extract double peaks to ensure the robustness and stability of bottom identification. Step 6: Based on the difference in the speed of light propagation at the water-air interface, establish a physical correction model for refractive index, correct the original water depth data for refraction, and obtain the true bottom elevation. Step 7: Perform median filtering smoothing on the water depth sequence for the entire orbit, and output a multi-dimensional parameter dataset including latitude and longitude, water surface height, water depth, original water depth and corrected water depth. Step 8: Integrate several orbital data, fuse and remove anomalies from all transit orbital data within a fixed time period, and generate continuous underwater topographic results based on interpolation; Step 9: Use several beam sonar measurements of underwater topography data to verify the obtained underwater topography, and compare the corrected seabed elevation with the measured data to evaluate the inversion accuracy.

2. The underwater topography inversion method based on altimeter satellites according to claim 1, characterized in that: In step 1, the time of the optical image data and the time of the altimeter satellite data are controlled within one month. The optical image data includes Landsat and Sentinel-2 optical satellite data.

3. The underwater topography inversion method based on altimeter satellites according to claim 1, characterized in that: In step 1, the normalized water index (NDWI) optical image is used to extract the water body vector boundary. The spectral bands used in the formula vary depending on the satellite data source. The formula for the normalized water index (NDWI) is as follows: In the formula This represents the green band, corresponding to the B2 band of Landsat 5 and Landsat 7, and the B3 band of Landsat 8, Landsat 9, and Sentinel-2. Represents the near-infrared band, corresponding to the B4 band of Landsat 5 and Landsat 7, the B5 band of Landsat 8 and Landsat 9, and the B8 band of Sentinel-2.

4. The underwater topography inversion method based on altimeter satellites according to claim 1, characterized in that: In step 2, the ICESat-2 ATL03 data is in h5 format. We select to extract strong beam photon information, combine it with the water body vector boundary to filter photon points within the water body, and output it in Shapefile format, which includes parameters such as latitude and longitude, ellipsoidal height, and signal confidence.

5. The underwater topography inversion method based on altimeter satellites according to claim 1, characterized in that: In step 4, during the photon signal processing stage, since the number of photons on the water surface is greater than the set value and is easy to identify, a random denoising strategy based on height distribution is adopted with 75 photons as the sliding window. The upper noise is reduced according to the height distribution, and the weak signal at the bottom is highlighted.

6. The underwater topography inversion method based on altimeter satellites according to claim 1, characterized in that: In step 6, the specific process of establishing the light propagation refractive index correction model is as follows: Taking into account the difference in propagation speed between air and water, the preliminary water depth results are corrected for refraction to obtain the true bottom elevation. The refraction correction model is based on the air refractive index. With the refractive index of water The ratio is corrected as follows: In the formula This is the corrected underwater elevation, where This refers to the water surface elevation. This represents the original underwater elevation.

7. The underwater topography inversion method based on altimeter satellites according to claim 1, characterized in that: In step 7, the underwater elevation sequence along the entire track is smoothed by median filtering to eliminate local outliers and form a continuous and stable underwater topographic profile. The dataset is then used for underwater topographic modeling and mapping.

8. The underwater topography inversion method based on altimeter satellites according to claim 1, characterized in that: In step 8, all transit orbit data from 2019 to 2025 are integrated, abnormal underwater elevations are removed, and a digital elevation model is generated based on interpolation methods.

9. The underwater topography inversion method based on altimeter satellites according to claim 1, characterized in that: In step 9, the accuracy of the inversion results is verified using underwater topographic data measured by multibeam sonar. During the verification process, attention should be paid to the correction of geoid anomalies between data. By comparing the corrected underwater elevation with the sonar measurement results, the inversion accuracy is evaluated, thereby ensuring the reliability of the method.

Images