An underwater terrain inversion method based on a spaceborne ICESat-2 photon satellite

Abstract

The application discloses a kind of underwater terrain inversion methods based on satellite ICESat-2 photon satellite, it is related to remote sensing water depth inversion technical field.Application ICESat-2 photon denoising and local robust estimation method is carried out inversion, data acquisition and pretreatment;Large-area photon denoising and signal extraction based on PAConv-PointNet++, PAConv-PointNet++ is for directly embedding PAConv into PointNet++;Water depth physical correction, including water depth refraction correction and water depth tide correction;Based on local robust estimation method, fine denoising and water depth trajectory fitting are carried out.The application does not need artificial preset clustering parameter, and is autonomously learned in a data-driven manner, in a complex shallow water environment of high turbidity, very shallow water, high noise ratio, signal photon and noise photon can still be accurately separated, and the denoising performance and environmental adaptability of sparse photon data are significantly improved.

Description

Technical Field

[0001] This invention relates to the field of remote sensing water depth inversion technology, and in particular to an underwater topography inversion method based on the spaceborne ICESat-2 photon satellite. Background Technology

[0002] Nearshore shallow waters are the core area of ​​land-sea interaction. Influenced by multiple factors such as waves, tides, and runoff, the topography changes dramatically, and the tidal channel system evolves frequently. Accurate measurement of water depth data and topographic dynamics in shallow waters is of great theoretical and practical significance for the study of coastal zone evolution mechanisms, marine resource management, navigation safety, and ecological environmental protection.

[0003] However, existing water depth measurement technologies face numerous challenges in complex shallow water environments. Traditional shipborne acoustic sounding, while achieving centimeter-level accuracy, suffers from high operating costs, limited coverage, and a high risk of grounding in shallow areas, making it difficult to reach remote or dangerous waters. Airborne lidar sounding compensates for the limitations of ships to some extent, but its technical threshold and cost per unit area are high, relying on dedicated software and professional personnel, making it difficult to achieve large-scale continuous sounding. The Advanced Terrain Laser Altimeter System (ATLAS) carried by the ICESat-2 satellite can emit 532nm wavelength green lasers, possessing the ability to penetrate water. Its ATL03 global geolocation photon data provides a new high-precision data source for nearshore optical shallow water water depth inversion. However, existing water depth denoising methods based on ICESat-2 photon data have significant technical defects: 1) Lasers undergo strong absorption and scattering effects during transmission through water, resulting in extremely sparse signal photon density in the bottom echo as water depth increases; daytime observation data is severely affected by solar background noise, with noise photons accounting for a much higher proportion than effective signal photons, leading to reduced water depth accuracy. 2) Traditional density clustering denoising methods such as DBSCAN and OPTICS require manual adjustment of parameters such as neighborhood radius and minimum number of points, resulting in poor parameter adaptability. They only perform well in clear water, and their denoising performance drops sharply in complex environments such as turbid or extremely shallow water. 3) Even after preliminary denoising, a large amount of noise and outliers remain in the underwater signal photons. Existing least squares and other fitting methods are highly sensitive to outliers and are easily affected by outlier interference, leading to distortion of the water depth fitting results. They lack robust outlier suppression and continuous water depth trajectory fitting methods.

[0004] In view of the above problems, existing methods are difficult to meet the requirements of accuracy, robustness and environmental adaptability of ICESat-2 photon data for water depth measurement in complex shallow water environments. Therefore, this invention is proposed. Summary of the Invention

[0005] To address the issues of low denoising accuracy, high degree of manual intervention, and poor scene adaptability of ICESat-2 photon data, this invention proposes a water depth retrieval method based on the spaceborne ICESat-2 photon satellite. This method is based on a two-step optical shallow water ICESat-2 photon denoising and local robust estimation water depth retrieval approach. It significantly improves the denoising performance and environmental adaptability of sparse photon data.

[0006] The objective of this invention is achieved through the following technical solutions:

[0007] An underwater topography inversion method based on the spaceborne ICESat-2 photon satellite, employing ICESat-2 photon denoising and local robustness estimation methods for inversion, includes the following steps:

[0008] S1: Data acquisition and preprocessing;

[0009] S2: Large-area photonic denoising and signal extraction based on the PAConv-PointNet++ model, where the PAConv-PointNet++ model is a direct embedding of PAConv into the PointNet++ model;

[0010] S3: Physical correction of water depth, including water depth refraction correction and water depth tidal correction;

[0011] S4: Refined denoising and water depth trajectory fitting based on local robust estimation method.

[0012] Preferably, step s1 includes:

[0013] S11: Collect ICESat-2 ATL03 global geolocation photon data of the study area, filter satellite orbit data that extend along the nearshore zone and include seabed photon reflections, and complete the batch download of raw data;

[0014] S12: Perform initial screening on the raw ATL03 data, extract ocean-type photons, and remove extreme outliers based on the effective depth measurement range of the optical shallow water area.

[0015] S13: Perform coordinate transformation on the filtered photon data, projecting the original latitude and longitude geographic coordinates to the UTM coordinate system;

[0016] S14: Segment the projected photon data;

[0017] S15: Collect tidal model data corresponding to the study area;

[0018] S16: Collect measured water depth verification data for the study area and unify coordinate and elevation benchmarks.

[0019] Preferably, step S2 includes:

[0020] S21: Based on the characteristics of the photons along the orbit exhibiting a double-layer structure consisting of sea surface photons and seabed photons in the elevation direction, a preliminary classification of the double-layer structure photons is made.

[0021] S22: Construct a PAConv-PointNet++ deep learning model to extract adaptive features of photon distribution and achieve large-area denoising of underwater photons.

[0022] Preferably, step S21 includes:

[0023] S211: By estimating the kernel density of photon elevation distribution along the track, the photon characteristics of the sea surface-seabed double-layer structure are quantitatively identified, the sea surface photon signal is marked, and the preliminary extraction of sea surface photons is completed.

[0024] S212: Based on the characteristics of continuous seabed topography and small elevation difference between adjacent seabed photons, calculate the elevation difference between adjacent seabed photons. Combining the elevation difference frequency histogram with the exponential decay model to screen for those that meet the requirements. Using the maximum elevation difference as a threshold, three rounds of processing are performed to progressively eliminate noisy photon points and mark the effective signal photons on the seabed. c y is a constant term, and y is the model fitting frequency corresponding to the elevation difference Δelev between adjacent photons;

[0025] S213: Based on the discrete characteristics of the elevation difference distribution of randomly distributed noise photons, photons that do not meet the elevation difference threshold are marked as noise photons, thus completing the preliminary division of signal photons and noise photons on the sea surface and seabed.

[0026] S214: Employing a lightweight strategy of fewer batches, multiple rounds, and modularization, it performs manual verification and correction on samples from complex scenarios. For the signal aliasing problem in extremely shallow water scenarios, it defines the physical boundary between photons on the sea surface and the seabed, and corrects erroneous divisions of aliasing regions. For the signal-to-noise boundary blurring problem in high-noise scenarios, it combines the optical characteristics of water bodies and the spatial distribution patterns of photons to distinguish weak signal photons from high-density noise clusters, avoiding the erroneous rejection of valid seabed signals.

[0027] Preferably, step S22 includes:

[0028] S221: In the PAConv-PointNet++ model, a sampling layer is used to sample the original input point cloud data, select the center point of the local neighborhood, and use the farthest point sampling method to ensure that the sampling center is evenly distributed in space.

[0029] S222: In the PAConv-PointNet++ model, the grouping layer is responsible for fusing the local regions of the sampling points. For the input sampling point matrix, each sampling point is used as the center of a circle. Combining the center of the circle with the region radius, a local region is constructed. Local feature points are sampled in each region to obtain a point set matrix divided into local regions.

[0030] S223: In the PAConv-PointNet++ model, PAConv is responsible for encoding the grouped local neighborhood features and outputting a discriminative local feature vector.

[0031] S224: Generate dynamic convolutional kernels that are adapted to the current point pair position using the three modules of PAConv in the PAConv-PointNet++ model: weight library, scoring network, and dynamic kernel.

[0032] S225: The feature stitching stage adopts a hierarchical and progressive implementation method, which integrates distance-weighted interpolation and cross-layer skip connection mechanism. The stitched features are input into the shared MLP unit for nonlinear transformation and feature extraction, and then propagated layer by layer to the original point cloud resolution. Finally, the classification probability of each photon point is output through the fully connected layer to complete the point-by-point segmentation task of signal photons and noise photons.

[0033] S226: Use the labeled training dataset to train and optimize the model. After training, input the underwater photon dataset into the model and output the classification probability of each photon point to complete the point-by-point segmentation of underwater signal photons and noise photons.

[0034] Preferably, water depth refraction correction includes:

[0035] Based on optical propagation characteristics and Snell's law, a refraction correction model for laser at the air-water interface is established. A quantitative relationship between the laser refraction angle and the incident angle is defined. Using trigonometric functions, a formula for calculating the correction amount for underwater photons in the vertical direction is derived. The formula is as follows:

[0036] ,

[0037] in, Q is Geometric deviation vector magnitude, In the formula, S is the uncorrected slant range, and R is the refraction-corrected slant range. , In the formula, and Let n1 and n2 represent the angle of incidence and the angle of refraction, respectively, where n1 is the refractive index of air and n2 is the refractive index of seawater. D This indicates the uncorrected water depth value.

[0038] Preferably, the correction amount is calculated by decomposing the angle to determine the offset direction. The formula for calculating the offset direction by decomposing the angle is as follows: , where γ is the total deviation angle, α is the vertical deviation angle, and β is the horizontal deviation angle.

[0039] Preferably, the water depth tidal correction includes:

[0040] The tidal level value corresponding to the observation time is obtained from the tidal model or tidal level observation data, and the water depth value under the unified datum is obtained by correction. The specific correction formula is as follows.

[0041] ,

[0042] in, To complete the tidal-corrected water depth value, The underwater photon elevation acquired by ICESat-2 at time t. This is the instantaneous tide level value. ,in, , and Let g represent the amplitude, angular frequency, and initial phase of the g-th tidal component, respectively, where g takes values ​​of 1, 2, 3...G, and G is the number of tidal components involved in the calculation.

[0043] Preferably, step S4 includes:

[0044] S41: Construct an iterative weighted model for local robust estimation. Using the assumption of equal-weighted independent observation, construct the error equation of the Gauss-Markov model, as shown in Equation 19. Combine the local robust estimation to construct the objective function, as shown in Equation 20. Introduce the IGG3 scheme based on the iteration of a posteriori unit weight error to construct the robust equivalent weight function, as shown in Equation 21. Continuously reduce the weight of outliers in the estimation process through the iterative weighting strategy.

[0045] 19,

[0046] in, V For the residual vector, A The coefficient matrix, L For the observation vector, P For the observation weight matrix, Let D(L) be the variance-covariance matrix of the observed vectors, representing the post-hoc unit weight variance.

[0047] 20,

[0048] in, For matrix A The i One portion, Let i be the i-th component of the observation vector L. The weight corresponding to the i-th observation. V i Let n be the i-th residual, and n be the total number of observations, i.e., the total number of photon points within the strip. For parameters The estimated value, The weight function in robust estimation, The operation refers to finding the minimum value;

[0049] twenty one,

[0050] in, For residual The standardized residual components, Equivalent rights are rights over the original observation rights. The new weights obtained after correction and Let be an empirical constant, and take . , The specific value can be adjusted according to the data size and noise characteristics;

[0051] S42: Combine the ICESat-2 system parameters, set the sliding window parameters, and use the median elevation value of the depth-sensing photons in each window as the initial parameter estimate;

[0052] S43: Perform local robust estimation iterative solution. Within each sliding window, update the equivalent weight matrix and the estimated model parameters iteratively until the iteration converges. The equivalent weights of the outlier observations are gradually reduced to close to zero, thus effectively suppressing outliers.

[0053] S44: Based on the iterative solution results, the discrete seabed signal photon points are transformed into spatially continuous and smooth water depth trajectories along the track, and the final water depth inversion results are output.

[0054] Compared with the prior art, the present invention has the following advantages:

[0055] First, leveraging the powerful spatial feature extraction capabilities of the PAConv-PointNet++ deep learning model, high-precision preliminary segmentation of signal photons and noise photons is achieved. Second, by integrating refraction correction and tidal correction, the refraction effect of laser at the air-water interface and the systematic errors caused by the dynamic tidal changes at the observation time are systematically eliminated, achieving a unified water depth benchmark. Finally, the second step of refined denoising and fitting is performed. Through the local robust estimation method, an iterative weighted model is constructed using the IGG3 weight function, which effectively suppresses and removes the outlier interference remaining in the bottom signal photons after the initial large-area denoising. This transforms discrete photon points into continuous and smooth water depth trajectories, solving the problems of discontinuous water depth inversion results and large outlier influence under sparse photon conditions, further improving the accuracy and stability of water depth inversion.

[0056] This invention employs a two-step photon denoising strategy combining kernel density estimation for preliminary classification with the PAConv-PointNet++ deep learning model. It eliminates the need for manually pre-setting clustering parameters and autonomously learns photon feature patterns through a data-driven approach. The trained PAConv-PointNet++ model achieves a recall of 95.59%, precision of 95.17%, F1 score of 95.38%, and IoU of 91.17%. Even in complex shallow water environments with high turbidity, extremely shallow water, and high noise levels, it can still achieve high-precision segmentation of signal and noise photons, significantly improving the denoising performance and environmental adaptability of sparse photon data. Attached Figure Description

[0057] Figure 1 This is an overall framework diagram of the present invention;

[0058] Figure 2 The images show the physical calibration diagram of photon depth using ICESat-2 ATL03 and the scatter plot of the depth extraction results from ATL03 compared with measured data. Among them, (a) is the physical calibration diagram of photon depth using the ATL03_20190530065739_09470307_007_gt3l strip, (b) is the scatter plot of the depth inversion results of this method compared with measured data, (c) is the physical calibration diagram of photon depth using the ATL03_20190228111804_09470207_007_gt3l strip, and (d) is the scatter plot of the depth inversion results of this method compared with measured data.

[0059] Figure 3The following are comparison figures of the denoising results of the PAConv-PointNet++ model and density clustering algorithms (OPTICS, DBSCAN): (a) is a schematic diagram of the original ICESat-2 data for a clear scene; (b) is a schematic diagram of the original ICESat-2 data for a very shallow scene; (c) is a schematic diagram of the original ICESat-2 data for a high-noise / low-transparency scene; (d) is the denoising result of the original ICESat-2 data for a clear scene using the OPTICS method; (e) is the denoising result of the original ICESat-2 data for a very shallow scene using the OPTICS method; (f) is the denoising result of the original ICESat-2 data for a high-noise / low-transparency scene using the OPTICS method; (g) is the denoising result of the original ICESat-2 data for a clear scene using the DBSCAN method; (h) is the denoising result of the original ICESat-2 data for a very shallow scene using the DBSCAN method; and (i) is the denoising result of the original ICESat-2 data for a high-noise / low-transparency scene using the DBSCAN method.

[0060] Figure 4 The following are the results of photon segmentation in complex environments for multiple scenes: (a) is a schematic diagram of the original ICESat-2 data for a clear scene, (b) is a schematic diagram of the labeled data for a clear scene, (c) is a water depth inversion result of the method in a clear scene; (d) is a schematic diagram of the original ICESat-2 data for an extremely shallow scene, (e) is a schematic diagram of the labeled data for an extremely shallow scene, (f) is a water depth inversion result of the method in an extremely shallow scene; (g) is a schematic diagram of the original ICESat-2 data for a high-noise / low-transparency scene, (h) is a schematic diagram of the labeled data for a high-noise / low-transparency scene, and (i) is a water depth inversion result of the method in a high-noise / low-transparency scene. Detailed Implementation

[0061] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0062] Implementation 1

[0063] This embodiment provides an underwater terrain inversion method based on the spaceborne ICESat-2 photon satellite. The inversion is performed using ICESat-2 photon denoising and local robustness estimation methods, and includes the following steps:

[0064] S1: Data Acquisition and Preprocessing

[0065] S11: Collect ICESat-2 ATL03 global geolocation photon data of the study area, filter satellite orbit data that extend along the nearshore zone and include seabed photon reflections, and complete the batch download of raw data;

[0066] S12: Perform initial screening on the raw ATL03 data, extract ocean-type photons, and remove photons from non-oceanic scenes such as land and sea ice; based on the effective depth sounding range of the optical shallow water area, retain photon data with elevation in the range of -50m to 50m, and remove extreme outliers.

[0067] S13: Perform coordinate transformation on the filtered photon data, projecting the original latitude and longitude geographic coordinates to the UTM coordinate system to eliminate the influence of geographic coordinate distortion on spatial feature extraction;

[0068] S14: The projected photon data is segmented. In this embodiment, the photon data of each orbit is divided into units of 8192 photon points, generating several subsets containing seabed photons (see...). Figure 4 This provides standardized input for subsequent model processing;

[0069] S15: Collect tidal model data corresponding to the study area, including TPXO and FES global / regional tidal model data, for subsequent tidal correction processing;

[0070] S16: Collect measured water depth verification data for the study area, unify coordinates and elevation benchmarks, and use them to verify the accuracy of subsequent inversion results.

[0071] S2: Large-area photon denoising and signal extraction based on the PAConv-PointNet++ model, where the PAConv-PointNet++ model is formed by directly embedding PAConv into the PointNet++ model; specifically, it includes the following steps:

[0072] S21: Based on the characteristics of the along-track photons exhibiting a double-layered structure consisting of sea surface photons and seabed photons in the elevation direction, a preliminary classification of the double-layered structure photons is performed; including the following steps:

[0073] S211: By performing kernel density estimation (KDE) on the elevation distribution of photons along the track (see Equation 1), the characteristics of the sea surface-seabed double-layer structure are quantitatively identified, sea surface signal photons are labeled, and the initial extraction of sea surface photons is completed.

[0074] (1),

[0075] in, x For the point whose elevation is to be estimated, Indicates the first The elevation of a photon This represents the total number of photons within the strip. This is the bandwidth parameter, used to control the smoothness of the density estimation. The Gaussian kernel function;

[0076] S212: Based on the characteristics of continuous seabed topography and small elevation difference between adjacent seabed photons, calculate the elevation difference between adjacent photons. (See Equation 2), combining the elevation difference frequency histogram and the exponential decay model (see Equation 3) to screen for those that meet the requirements. The maximum elevation difference under the condition is used as a threshold, and three rounds are performed to gradually eliminate noisy photon points and mark the effective signal photons on the seabed;

[0077] (2),

[0078] (3),

[0079] in It is the height of each photon. a It is the amplitude of the fitted decay model. b The attenuation rate, c For constant terms;

[0080] S213: Based on the discrete characteristics of the elevation difference distribution of randomly distributed noise photons, photons that do not meet the elevation difference threshold are marked as noise photons, thus completing the preliminary division of signal photons and noise photons on the sea surface and seabed.

[0081] S214: Employing a lightweight strategy of fewer batches, multiple rounds, and modularization, it performs manual verification and correction on samples from complex scenarios. For the signal aliasing problem in extremely shallow water scenarios, it defines the physical boundary between sea surface and seabed photons based on prior knowledge of continuous changes in seabed topography, correcting erroneous divisions of aliasing regions. For the signal-to-noise boundary blurring problem in high-noise scenarios, it distinguishes weak signal photons from high-density noise clusters by combining the optical characteristics of water bodies and the spatial distribution patterns of photons, avoiding the erroneous rejection of valid seabed signals.

[0082] S22: Construct the PAConv-PointNet++ deep learning model to adaptively extract features from photon distribution and achieve large-area denoising of underwater photons, including the following steps:

[0083] S221: In the PAConv-PointNet++ model, the sampling layer samples the original input point cloud data and selects the center point of the local neighborhood. Considering the sparse and unevenly distributed photons in ICESat-2, this embodiment uses the farthest point sampling (FPS) method to ensure that the sampling centers are uniformly distributed in space, avoiding sampling deviations caused by signal photon aggregation and noise photon dispersion in complex shallow water areas. First, a starting point is selected, and the distance between the starting point and other points is calculated. The point with the largest distance is added to the downsampling point set, and then the distance is updated (see Equation 4). Then, the update is iterated downwards until the target number of sampling points is obtained. (See Formula 5) By maximizing the minimum distance criterion, the maximum coverage of sampling points in three-dimensional space is guaranteed;

[0084] (4)

[0085] (5)

[0086] in For the first m Layer t One sampling center, Let x be the input point set of the (m-1)th layer. Any point in, It is the L2 norm (Euclidean distance). , Let m be the set of sampling centers. Let be the number of sampling centers in the m-th layer. Let j be the index of the selected sampling centers that have been selected during the sampling process in the m-th layer. Solving for the operator at the location of the maximum value Nearest neighbor distance filtering operator.

[0087] S222: In the PAConv-PointNet++ model, the grouping layer is responsible for fusing local regions of the sampling points. For the input sampling point matrix, each sampling point is used as the center of a circle. Combining this center with a specific region radius, a local region is constructed. Within each region, local feature points are sampled to obtain a point set matrix divided into local regions. This embodiment uses a sphere query combined with a multi-scale grouping (MSG) strategy to adapt to the multi-scale structural features of photons in complex shallow water. For each sampling center... and preset multi-scale radius Building a neighborhood (See Equation 6), where the multi-scale radii in this embodiment are respectively , Then, the relative coordinates of each photon point in the neighborhood are normalized (see Equation 7) to eliminate the interference of absolute coordinate offset on feature extraction. This operation ensures that feature extraction focuses on local geometry rather than global position, enhancing the translation invariance of the model;

[0088] (6)

[0089] (7)

[0090] in, For the input point set Any point in the sampling center , rLet be the radius of the spherical neighborhood.

[0091] S223: In the PAConv-PointNet++ model, the PointNet layer is the core of feature extraction at each level, responsible for encoding the grouped local neighborhood features and outputting discriminative local feature vectors. Its design strictly adheres to the permutation invariance of point cloud data. However, the standard PointNet layer uses a multilayer perceptron (MLP) to perform the same feature transformation on all points. This position-independent processing method cannot adapt to the complex and varied geometric structures in point clouds. In complex shallow water environments, the undulating seabed topography and the random distribution of noise photons require the model to have position awareness capabilities. Therefore, this study introduces PAConv. PAConv constructs position-adaptive convolutional kernels through a dynamic assembly strategy, enabling the network to flexibly learn and respond to the unique geometric patterns of local regions. PAConv is then directly embedded into the PointNet++ model.

[0092] S224: The PAConv-PointNet++ model utilizes three modules—weight library, scoring network, and dynamic kernel generation—to generate a dynamic convolutional kernel adapted to the current point pair's position. In this embodiment, the weight library is defined as... Each element The weights are learnable weight matrices, where M is the number of weight matrices. The scoring network learns a non-linear mapping between the positional relationships of the point cloud and the matrices in the weight library, assigning adaptive coefficients to each weight matrix. The input to the scoring network (ScoreNet) is a positional relationship encoding vector containing relative coordinates. The absolute coordinates of the center point and its neighboring points and Euclidean distance The input format, which integrates full-dimensional location information in three-dimensional space, can optimally utilize spatial correlation information. ScoreNet outputs a normalized coefficient vector to quantify the contribution of each fundamental matrix in kernel assembly. Dynamic kernel generation enables the kernel function to adapt to local spatial features of the point cloud, flexibly modeling irregular geometric structures through a data-driven approach. The coefficients output by ScoreNet are used to perform a weighted summation of the fundamental matrices in the weight library (see Equation 8).

[0093] (8),

[0094] in Corresponding weight matrix The adaptive coefficient.

[0095] S225: The feature stitching stage adopts a hierarchical and progressive implementation approach, integrating distance-weighted interpolation and cross-layer skip connection mechanisms. In this embodiment, the feature value f of the sampling point is interpolated to the coordinates of all points in the previous layer based on the k-nearest neighbor inverse distance-weighted interpolation method (see Equation 9). To compensate for the details that may be lost during the interpolation process, the interpolated features are stitched together with the skip connection features from the set abstraction layer, fully integrating feature information from different levels. The stitched features are input into a shared MLP unit for nonlinear transformation and feature extraction, propagating layer by layer to the original point cloud resolution. Finally, the classification probability of each photon point is output through a fully connected layer, completing the point-by-point segmentation task of signal photons and noise photons.

[0096] (9),

[0097] in j is the feature channel index. Inverse distance weights Indicates the target point With the Neighboring points The distance between them Let x be the j-th dimension of the interpolated eigenvalue. For the first i The nearest neighbor point at the th j Known feature values ​​on each feature channel;

[0098] S226: Use the labeled training dataset to train and optimize the model. After training, input the underwater photon dataset into the model and output the classification probability of each photon point to complete the point-by-point segmentation of underwater signal photons and noise photons.

[0099] S3: Physical correction of water depth, including the following steps:

[0100] S31: Water depth refraction correction, including the following steps:

[0101] S311: When an ATLAS laser pulse transitions from air to seawater, a refraction effect occurs, thus requiring refraction correction for seabed photons. Based on optical propagation characteristics and Snell's law (see Equation 10), a refraction correction model for lasers at the air-water interface is established. The quantitative relationship between the laser refraction angle (see Equation 10) and the incident angle (see Equation 11) is defined, and the calculation formulas for the horizontal and vertical correction amounts of underwater photons are derived by combining trigonometric functions.

[0102] (10),

[0103] (11),

[0104] in and These represent the angle of incidence and the angle of refraction, respectively. ref_elev n1 is the incident elevation angle of the beam within the segment, n2 is the refractive index of air, and n3 is the refractive index of seawater.

[0105] S312: Based on the optical propagation characteristics of water, after the laser is refracted at the air-water interface, the uncorrected slant distance is S (see Equation 12) and the refracted slant distance is R (see Equation 13).

[0106] (12),

[0107] (13),

[0108] in, D This indicates the uncorrected water depth value.

[0109] S313: Geometric deviation vector magnitude of the laser path determined by the law of cosines Q (See Equation 14);

[0110] (14),

[0111] S314: Calculate the offset direction by angular decomposition (see Equation 15).

[0112] (15), where γ is the total deviation angle, α is the vertical deviation angle, and β is the horizontal deviation angle.

[0113] S315: Since the ICESat-2 satellite observes the Earth almost vertically, the incident elevation angle is very small. At a water depth of 30m, the horizontal correction is 9cm. The error in the horizontal direction is very small and can be ignored. In this embodiment, only the vertical refraction correction is calculated (see Equation 16).

[0114] (16),

[0115] S32: Tidal correction for water depth, including the following steps:

[0116] S321: The photon elevation acquired by the ICESat-2 satellite is essentially the instantaneous water level height relative to a reference ellipsoid or geoid at the observation time. Influenced by astronomical tides, meteorological conditions, and regional hydrodynamic processes, sea level varies significantly in time and space; therefore, tidal effect correction is performed. The tidal level value corresponding to the observation time is obtained through a tidal model or tidal observation data and removed from the observed elevation to obtain the water depth value under a unified reference surface. The specific correction formula is shown in Equation 17.

[0117] (17),

[0118] in, To complete the tidal-corrected water depth value, The underwater photon elevation acquired by ICESat-2 at time t.

[0119] S322: Based on the geographical location and observation time of the underwater photon, the instantaneous tide level value at the corresponding location and time is obtained by interpolation from the TPXO / FES tidal model. (See Equation 18) This type of model, based on harmonic analysis, represents tides as a linear superposition of several major tidal components, in the form of:

[0120] (18),

[0121] in, , and Let g represent the amplitude, angular frequency, and initial phase of the g-th tidal component, respectively, where g takes values ​​of 1, 2, 3, ..., G, and G is the number of tidal components involved in the calculation.

[0122] S4: Refined denoising and water depth trajectory fitting based on local robust estimation, including the following steps:

[0123] S41: To address the outliers and data sparsity issues in ICESat-2 laser altimetry satellite data, an iterative weighted model for local robust estimation is constructed. Using the assumption of equal-weighted independent observations, the error equation of the Gauss-Markov model is constructed (see Equation 19). Combined with local robust estimation, an objective function is constructed (see Equation 20). The IGG3 scheme based on a posteriori unit weight error iteration is introduced to construct a robust equivalent weight function (see Equation 21). The weight of outliers in the estimation process is continuously reduced through an iterative weighting strategy.

[0124] (19),

[0125] in, V For the residual vector, A The coefficient matrix, L For the observation vector, P For the observation weight matrix, Let D(L) be the variance-covariance matrix of the observed vectors, representing the post-hoc unit weight variance.

[0126] (20),

[0127] in, For matrix A The i One portion, Let i be the i-th component of the observation vector L. The weight corresponding to the i-th observation. V i Let n be the i-th residual, and n be the total number of observations, i.e., the total number of photon points within the strip. For parameters The estimated value, Let be the weight function in robust estimation, and i be the index of the independent photon point obtained by local robust estimation, taking values ​​from 1, 2, 3, ..., n. The operation refers to finding the minimum value;

[0128] (twenty one),

[0129] in, For residual The standardized residual components, Equivalent rights are rights over the original observation rights. The new weights obtained after correction and Let be an empirical constant, and take . , The specific value can be adjusted according to the data size and noise characteristics;

[0130] S42: Based on the ICESat-2 system parameters, set the sliding window parameters, and use the median elevation value of the depth-sensing photons within each window as the initial parameter estimate. Since the spot diameter of ICESat-2 is approximately 12m and the interval between adjacent laser pulses along the track is approximately 0.7m, this embodiment sets the length of the sliding window to 12m and the sliding step size to 0.7m.

[0131] S43: Perform local robust estimation iterative solution. Within each sliding window, update the equivalent weight matrix (see Equation 22) and the estimated model parameters iteratively until the iteration converges. The equivalent weights of the corresponding outlier observations are gradually reduced to close to zero, thereby effectively suppressing outliers.

[0132] (twenty two),

[0133] in, For the first n Equivalent weighting factor for each observation, It is an equivalent weight matrix.

[0134] S44: Based on the iterative solution results, the discrete seabed signal photon points are transformed into spatially continuous and smooth water depth trajectories along the track, and the final water depth inversion results are output.

[0135] Verification Example 1

[0136] To systematically verify the effectiveness and generalization ability of the PAConv-PointNet++ model in the seabed topography segmentation task, accuracy, precision, recall, F1 score, and intersection-over-union ratio (IoU) were used as evaluation indicators (see Table 1) to comprehensively evaluate the accuracy of the method in seabed signal extraction.

[0137] Table 1 Statistical Indicators

[0138]

[0139] Where TP is the number of true positives, FP is the number of false positives, and FN is the number of false negatives. This is the water depth value of the i-th photon point after ICESat-2 data correction. This is the in-situ depth sounding value corresponding to the i-th photon point, where i takes values ​​of 1, 2, 3…n, and n is the total number of photon points within the strip. Precision measures the proportion of correctly identified predictions; accuracy represents the proportion of samples predicted as belonging to a certain class that actually belong to that class, reflecting the model's reliability in its output; recall measures the proportion of samples that are successfully identified out of all samples that actually belong to that class, reflecting the model's coverage capability; the F1 score comprehensively balances precision and recall, using a harmonic average to unify the evaluation metrics, which is particularly important in applications with imbalanced multi-class distributions; IoU measures the degree of matching between the predicted results and the true values, simultaneously reflecting the model's accuracy in identifying target classes / target regions and its boundary fitting ability.

[0140] This invention employs a two-step photon denoising strategy combining kernel density estimation for preliminary classification with the PAConv-PointNet++ deep learning model. It eliminates the need for manually pre-setting clustering parameters and autonomously learns photon feature patterns through a data-driven approach. The trained PAConv-PointNet++ model achieves a recall of 95.59%, precision of 95.17%, F1 score of 95.38%, and IoU of 91.17%. Even in complex shallow water environments with high turbidity, extremely shallow water, and high noise levels, it can still achieve high-precision segmentation of signal and noise photons, significantly improving the denoising performance and environmental adaptability of sparse photon data.

[0141] Inversion accuracy verification involves matching and comparing the final water depth inversion results with the measured water depth data, and calculating the mean square error (MSE). RMSE Mean absolute error ( MAE Precision evaluation indicators such as )

[0142] Conclusion: The depth values ​​extracted by the PAConv-PointNet++ deep learning model show significant consistency with the measured values ​​from LiDAR (see...). Figure 2 Coefficient of determination R2 The accuracy reached above 0.95, and the mean absolute error (MAE) was below 0.3m, indicating that the model output was close to unbiased in relation to the actual depth and had a good linear correlation. The accuracy evaluation results show that the PAConv-PointNet++ model has good accuracy and stability in the seabed photon segmentation task.

[0143] A comprehensive comparative analysis of the denoising results of the PAConv-PointNet++ model and density clustering algorithms (OPTICS, DBSCAN) is conducted to verify the denoising performance of the proposed method.

[0144] Conclusion: Through quantitative comparison with traditional OPTICS and DBSCAN density clustering algorithms, it was found that the denoising effect of traditional algorithms is affected by parameter settings, and they cannot simultaneously adapt to complex scenarios such as signal aliasing in extremely shallow water and high-noise, weak signals, making it difficult to achieve high-precision batch denoising across all scenarios (see...). Figure 3 The PAConv-PointNet++ model significantly outperforms traditional density clustering algorithms in three typical shallow sea scenarios: extremely shallow water, clear water, and high-noise / low-transparency water, with precision, recall, and F1 score all consistently exceeding 96%.

[0145] Sensitivity analysis was conducted to examine the denoising effect and inversion accuracy of the proposed method under different photon sparsity and noise environments, thus verifying the method's environmental adaptability.

[0146] Conclusion: By employing a locally robust estimation method and constructing an iterative weighted model using the IGG3 weight function, the residual outlier interference in the underwater signal photons is effectively suppressed. This transforms discrete photon points into continuous and smooth water depth trajectories, solving the problems of discontinuous water depth inversion results and significant outlier influence under sparse photon conditions. This further improves the accuracy and stability of water depth inversion (see...). Figure 4 ).

[0147] In summary, this invention constructs a two-step method for water depth inversion in optical shallow water areas using ICESat-2 photon denoising and local robust estimation. The first step combines kernel density estimation with the PAConv-PointNet++ deep learning model to overcome the shortcomings of traditional methods, such as high manual intervention and poor scene adaptability. This achieves large-area, high-precision segmentation of seabed signal photons from background noise in complex aquatic environments, followed by a unified benchmark through refraction and tidal correction. The second step uses a local robust estimation method and the IGG3 weighting function to refine the denoising of outliers remaining after the large-area denoising in the first step, effectively suppressing local outlier interference. Finally, the discrete points are fitted into a continuous and smooth water depth trajectory. This two-step process has a high degree of standardization, perfectly adapting to the complex environment of nearshore optical shallow water areas. It solves the problems of low water depth inversion accuracy and poor continuity under extremely sparse photon conditions, providing high-precision, large-scale water depth data support for dynamic topographic monitoring in nearshore waters.

[0148] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for underwater terrain inversion based on a spaceborne ICESat-2 photon satellite, characterized by, The inversion is performed using the ICESat-2 photonic denoising and local robustness estimation method, including the following steps: S1: Data acquisition and preprocessing; S2: Large-area photonic denoising and signal extraction based on the PAConv-PointNet++ model, where the PAConv-PointNet++ model is a direct embedding of PAConv into the PointNet++ model; including: S21: Based on the characteristic of photons exhibiting a two-layer structure consisting of surface photons and seabed photons along the elevation direction, a preliminary classification of two-layer structure photons is performed; including: S211: By estimating the kernel density of photon elevation distribution along the track, the photon characteristics of the sea surface-seabed double-layer structure are quantitatively identified, the sea surface photon signal is marked, and the preliminary extraction of sea surface photons is completed. S212: Based on the characteristics of continuous seabed topography and small elevation difference between adjacent seabed photons, calculate the elevation difference between adjacent seabed photons. Combining the elevation difference frequency histogram with the exponential decay model to screen for those that meet the requirements. Using the maximum elevation difference as a threshold, three rounds of processing are performed to progressively eliminate noisy photon points and mark the effective signal photons on the seabed. c y is a constant term, and y is the model fitting frequency corresponding to the elevation difference Δelev between adjacent photons; S213: Based on the discrete characteristics of the elevation difference distribution of randomly distributed noise photons, photons that do not meet the elevation difference threshold are marked as noise photons, thus completing the preliminary division of signal photons and noise photons on the sea surface and seabed. S214: Employing a lightweight strategy of fewer batches, multiple rounds, and modularization, it manually verifies and corrects samples from complex scenarios. For signal aliasing issues in extremely shallow water scenarios, it defines the physical boundaries between photons on the sea surface and the seabed, correcting erroneous divisions of aliasing regions. For the problem of blurred signal-to-noise boundaries in high-noise scenarios, it combines the optical characteristics of water bodies and the spatial distribution patterns of photons to distinguish weak signal photons from high-density noise clusters, avoiding the erroneous rejection of valid seabed signals. S22: Construct a PAConv-PointNet++ deep learning model to extract adaptive features of photon distribution and achieve large-area denoising of underwater photons; S3: Physical correction of water depth, including water depth refraction correction and water depth tidal correction; S4: Refined denoising and water depth trajectory fitting based on local robustness estimation method; including: S41: Construct an iterative weighted model for local robust estimation. Using the assumption of equal-weighted independent observation, construct the error equation of the Gauss-Markov model. Combine the local robust estimation to construct the objective function. Introduce the IGG3 scheme based on the iteration of a posteriori unit weight error to construct a robust equivalent weight function. Continuously reduce the weight of outliers in the estimation process through the iterative weighting strategy. The error equation is: , in, V For the residual vector, A The coefficient matrix, L For the observation vector, P For the observation weight matrix, Let D(L) be the variance-covariance matrix of the observed vectors, representing the post-hoc unit weight variance. The objective function is: , in, For matrix A The i One portion, Let i be the i-th component of the observation vector L. The weight corresponding to the i-th observation. V i Let n be the i-th residual, and n be the total number of observations, i.e., the total number of photon points within the strip. For parameters The estimated value, The weight function in robust estimation, The operation refers to finding the minimum value; The robust equivalent weight function is: , in, For residual The standardized residual components, Equivalent rights are rights over the original observation rights. The new weights obtained after correction and Let be an empirical constant, and take . , The specific value should be adjusted based on the data size and noise characteristics. S42: Combine the ICESat-2 system parameters, set the sliding window parameters, and use the median elevation value of the depth-sensing photons in each window as the initial parameter estimate; S43: Perform local robust estimation iterative solution. Within each sliding window, update the equivalent weight matrix and the estimated model parameters iteratively until the iteration converges. The equivalent weights of the outlier observations are gradually reduced to close to zero, thus effectively suppressing outliers. S44: Based on the iterative solution results, the discrete seabed signal photon points are transformed into spatially continuous and smooth water depth trajectories along the track, and the final water depth inversion results are output.

2. The underwater topography inversion method according to claim 1, characterized in that, Step s1 includes: S11: Collect ICESat-2 ATL03 global geolocation photon data of the study area, filter satellite orbit data that extend along the nearshore zone and include seabed photon reflections, and complete the batch download of raw data; S12: Perform initial screening on the raw ATL03 data, extract ocean-type photons, and remove extreme outliers based on the effective depth measurement range of the optical shallow water area. S13: Perform coordinate transformation on the filtered photon data, projecting the original latitude and longitude geographic coordinates to the UTM coordinate system; S14: Segment the projected photon data; S15: Collect tidal model data corresponding to the study area; S16: Collect measured water depth verification data for the study area and unify coordinate and elevation benchmarks.

3. The underwater topography inversion method according to claim 1, characterized in that, Step S22 includes: S221: In the PAConv-PointNet++ model, a sampling layer is used to sample the original input point cloud data, select the center point of the local neighborhood, and use the farthest point sampling method to ensure that the sampling center is evenly distributed in space. S222: In the PAConv-PointNet++ model, the grouping layer is responsible for fusing the local regions of the sampling points. For the input sampling point matrix, each sampling point is used as the center of a circle. Combining the center of the circle with the region radius, a local region is constructed. Local feature points are sampled in each region to obtain a point set matrix divided into local regions. S223: In the PAConv-PointNet++ model, PAConv is responsible for encoding the grouped local neighborhood features and outputting a discriminative local feature vector. S224: Generate dynamic convolutional kernels that are adapted to the current point pair position using the three modules of PAConv in the PAConv-PointNet++ model: weight library, scoring network, and dynamic kernel. S225: The feature stitching stage adopts a hierarchical and progressive implementation method, which integrates distance-weighted interpolation and cross-layer skip connection mechanism. The stitched features are input into the shared MLP unit for nonlinear transformation and feature extraction, and then propagated layer by layer to the original point cloud resolution. Finally, the classification probability of each photon point is output through the fully connected layer to complete the point-by-point segmentation task of signal photons and noise photons. S226: Use the labeled training dataset to train and optimize the model. After training, input the underwater photon dataset into the model and output the classification probability of each photon point to complete the point-by-point segmentation of underwater signal photons and noise photons.

4. The underwater topography inversion method according to claim 1, characterized in that, Water depth refraction correction includes: Based on optical propagation characteristics and Snell's law, a refraction correction model for laser at the air-water interface is established. A quantitative relationship between the laser refraction angle and the incident angle is defined. Using trigonometric functions, a formula for calculating the correction amount for underwater photons in the vertical direction is derived. The formula is as follows: ， in, Q is Geometric deviation vector magnitude, , , In the formula, S is the uncorrected slant range, and R is the refraction-corrected slant range. and Let n1 and n2 represent the angle of incidence and the angle of refraction, respectively, where n1 is the refractive index of air and n2 is the refractive index of seawater. D This indicates the uncorrected water depth value.

5. The underwater topography inversion method according to claim 4, characterized in that, The formula for calculating the correction amount is obtained by calculating the offset direction through angle decomposition. The formula for calculating the offset direction through angle decomposition is as follows: , where γ is the total deviation angle, α is the vertical deviation angle, and β is the horizontal deviation angle.

6. The underwater topography inversion method according to claim 1, characterized in that, Tidal depth corrections include: The tidal level value corresponding to the observation time is obtained from the tidal model or tidal level observation data. The water depth value under the unified datum is obtained by correction. The specific correction formula is as follows: ， in, To complete the tidal-corrected water depth value, The underwater photon elevation acquired by ICESat-2 at time t. This is the instantaneous tide level value. ,in, , and Let g represent the amplitude, angular frequency, and initial phase of the g-th tidal component, respectively, where g takes values ​​of 1, 2, 3, ..., G, and G is the number of tidal components involved in the calculation.

Images