Multi-temporal ICESat-2 turbid coastal water depth retrieval method and apparatus

Abstract

The present application discloses a turbid shallow sea depth inversion method of multi-time phase ICESat-2, comprising the following steps: extracting sea surface photon clusters from original photon point cloud data, denoising, segmenting and reconstructing sea surface height program sequence, extracting discrete water depth values of different water depth observation points in the nearshore shallow water research area by using an adaptive window, and based on the discrete water depth values of the water depth observation points, carrying out water depth inversion and multi-time phase data fusion on each target grid point covering the nearshore shallow water research area. The present application also proposes a turbid shallow sea depth inversion device of multi-time phase ICESat-2 for realizing the method. The present application overcomes the defect of inaccurate extraction of fixed window spectrum, and proposes a multi-time phase sounding point fusion method of segmented double weight distance from shore, generates a continuous underwater terrain digital elevation model DEM through multi-time phase data fusion, and realizes high-precision water depth inversion in turbid water area. The present application provides an important technical means for marine resource development and coastal zone management.

Description

Technical Field

[0001] This invention belongs to the field of spaceborne lidar marine remote sensing technology, specifically involving a method for inverting the depth of turbid shallow waters using multi-temporal ICESat-2, and also involving a device for inverting the depth of turbid shallow waters using multi-temporal ICESat-2, which is suitable for underwater topographic exploration of coastal zones. Background Technology

[0002] Coastal topographic data, including shallow water underwater topographic data, is traditionally acquired primarily through shipborne sonar bathymetry (such as single-beam / multi-beam echo sounders) and airborne lidar. While these methods can acquire high-precision data, they are costly, have limited coverage, and are difficult to implement for large-scale continuous monitoring in remote areas. With the development of satellite remote sensing technology, water depth retrieval methods based on optical imagery, synthetic aperture radar (SAR), and lidar have become cost-effective alternatives, but still have many limitations.

[0003] Optical satellite data is used to invert water depth using a radiative transfer model, based on the exponential decay of light in water. This method works well in clear waters, but in turbid waters, suspended sediments and organic matter significantly enhance light attenuation, preventing signals from penetrating the water and rendering the optical method almost ineffective. Water depth inversion methods based on synthetic aperture radar (SAR) backscattering intensity utilize the modulation of sea surface roughness by seabed topography to invert water depth through radar backscattering differences. This method is unaffected by water quality but relies on real-time ocean dynamic parameters such as tides and wind fields, and requires initial water depth data as prior knowledge. Furthermore, the visibility of topography in SAR images is limited by sea state conditions; weak signals under low wind speeds or complex flow fields lead to inversion failure. Satellite bathymetry based on wave theory inverts water depth by analyzing the deformation of surface waves in shallow water (wavelength shortening, wave velocity reduction). Its advantage lies in its independence from water transparency and the absence of prior topographic data.

[0004] However, existing water depth inversion methods based on ICESat-2 wave characteristics mostly employ fixed-window Fourier transforms, which cannot adapt to the dynamic changes in wave spectra. Especially in nearshore shallow waters, wave wavelengths shorten with decreasing water depth, and fixed windows may lead to incomplete or distorted spectral information extraction. Furthermore, ICESat-2 along-track data has sparse spatial coverage, making it difficult to generate a continuous underwater topographic digital elevation model (DEM) from single transit data, necessitating an effective multi-temporal data fusion method. To address these issues, this invention proposes a multi-temporal ICESat-2 water depth inversion method for turbid shallow waters. By optimizing wave spectrum extraction and multi-temporal data fusion, it significantly improves the accuracy and reliability of water depth inversion in turbid shallow waters. Summary of the Invention

[0005] The purpose of this invention is to address the problems of inaccurate fixed-window spectrum extraction and difficulties in multi-temporal data fusion in existing technologies. This invention proposes a multi-temporal ICESat-2 method for deriving turbid shallow water depth, and also presents a multi-temporal ICESat-2 turbid shallow water depth deriving device. This invention optimizes wave spectrum information extraction through adaptive window Fast Fourier Transform (FFT) and combines it with a dual-weighted fusion algorithm based on offshore distance segments to achieve the fusion of multi-temporal water depth results, generating a continuous underwater topographic digital elevation model (DEM).

[0006] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0007] The multi-temporal ICESat-2 method for inverting the depth of turbid shallow waters includes the following steps:

[0008] Step 1: Acquire multi-temporal ICESat-2 ATL03 level sea surface wave data covering the nearshore shallow water research area and the corresponding deep water area; extract raw photon point cloud data from the multi-temporal ICESat-2 ATL03 level sea surface wave data, and segment the photons in the raw photon point cloud data by the track direction and the vertical direction; for each track direction segment: based on the distribution of the number of photons in the vertical direction segment, select the corresponding sea surface photon clusters.

[0009] Step 2: Denoise the photon clusters on the sea surface.

[0010] Step 3: The space between adjacent sea surface photons in the denoised sea surface photon cluster that is greater than a preset gap threshold is used as the gap. Based on the gap, the sea surface photons corresponding to each beam are divided into multiple continuous segments. For each continuous segment, based on the along-track distance and elevation of the corresponding sea surface photon, smooth B-spline fitting and resampling are performed independently to generate a reconstructed sea surface elevation sequence. .

[0011] Step 4: Based on the reconstructed sea surface elevation sequence Using Fast Fourier Transform (FFT), the wave wavelengths corresponding to the waves with the largest amplitudes are extracted and corrected in the nearshore shallow water study area and the deep water area, respectively. Based on the corrected wave wavelengths, discrete water depth values ​​of different water depth observation points in the nearshore shallow water study area are extracted.

[0012] Step 5: Create a planar grid covering the nearshore shallow water study area, calculate and segment the offshore distances of all water depth observation points and target grid points in the planar grid; calculate the interpolation weights of the neighboring water depth observation points centered on each target grid point, and perform water depth interpolation calculations on all target grid points based on the interpolation weights of the neighboring water depth observation points and the discrete water depth values ​​to generate an underwater topographic digital elevation model (DEM).

[0013] The screening of sea surface photon clusters as described above specifically includes the following steps:

[0014] Step 1.1: Extract photon elevations from the original photon point cloud data, retaining photons whose elevation values ​​are within the preset elevation range, and use them as initial photons.

[0015] Step 1.2: For each beam in the original photon point cloud data, based on the beam trajectory, the photons along the track direction in the initial selection photons are segmented according to the preset track segment length, and each segment is recorded as a track direction segment.

[0016] Step 1.3: For each segment along the track direction: According to the preset vertical box height, the photons in the segment along the track direction are divided into vertical boxes. Each vertical box is called a vertical box.

[0017] Step 1.4: For each segment along the track direction, adaptively set the corresponding elevation extraction range of the sea surface photon cluster. ; This represents the height of the upper boundary of the vertical bin with the highest photon count. The lower boundary height of the vertical bin with the highest photon count. This indicates the height difference of the lower buffer zone. Indicates the height difference of the upper buffer zone:

[0018] .

[0019] .

[0020] in, This refers to the vertical compartment height; This represents the number of photons in the vertical bin with the largest number of photons. This represents the number of photons in the first vertical bin immediately below the vertical bin with the largest number of photons. This represents the number of photons in the first vertical box immediately above the vertical box with the largest number of photons.

[0021] The elevation values ​​of each segment along the track direction fall within the corresponding intervals. Photons within the ocean are considered as surface photons, and the collection of surface photons is called a surface photon cluster.

[0022] As described above, step 2 sequentially employs a density-based spatial clustering denoising algorithm and a two-dimensional sliding window midpoint filter to denoise the sea surface photons in the sea surface photon cluster.

[0023] As described above, step 3 specifically includes the following steps:

[0024] Step 3.1: For each beam, based on the beam trajectory, calculate the along-track spacing between adjacent sea surface photons in the denoised sea surface photon cluster, and use the space where the along-track spacing is greater than the preset gap threshold as the gap.

[0025] Step 3.2: Using the gaps identified in Step 3.1 as boundaries, divide the sea surface photons in the denoised sea surface photon cluster of Step 2 into multiple continuous segments; for each continuous segment, perform smooth B-spline fitting independently based on the orbital distance and elevation of the corresponding sea surface photon to generate the corresponding B-spline fitting curve function.

[0026] Step 3.3: For each beam: For each continuous segment corresponding to the beam, based on the number of evaluation points... , For the first The number of photons within a continuous segment is used to resample the corresponding B-spline fitting curve, resulting in the reconstructed sub-sea surface height sequence. All reconstructed sub-sea surface height sequences are then concatenated to obtain the reconstructed sea surface height sequence. .

[0027] As described above, step 4 specifically includes the following steps:

[0028] Step 4.1: Follow the reconstructed sea surface elevation sequence The corresponding resampling profile points are located in nearshore shallow water or deep water areas, which will reconstruct the sea surface elevation sequence. Divided into reconstructed nearshore sea surface elevation series and the reconstruction of deep-water sea surface elevation sequence .

[0029] Step 4.2: Serialize the reconstructed deep-water sea surface elevation for each beam. Perform a Fast Fourier Transform (FFT) to extract the wave wavelength corresponding to the wave with the maximum amplitude, and denote it as the main wave wavelength in the deep water area. .

[0030] Based on the prior initial water depth, an adaptive window is initialized; for each beam: the adaptive window slides along the beam trajectory of the nearshore shallow water study area based on a preset sliding step size; the reconstructed nearshore sea surface elevation is then programmed. The size of the adaptive window is iteratively adjusted using Fast Fourier Transform (FFT) and a preset wave wavelength range for the subsequences corresponding to each location of the adaptive window, and the main wave wavelength within the adjusted adaptive window range is extracted. The center of the adaptive window is the water depth observation point; This is the serial number of the water depth observation point.

[0031] Step 4.3: Use the observation phase difference between strong and weak beams to correct the angle between the wave propagation direction and the beam trajectory direction, thereby correcting the corresponding dominant wave wavelength. The dominant wave wavelength includes the dominant wave wavelength in deep water. and the main wavelength of the window The corrected deep-water wave wavelengths were obtained accordingly. and the corrected window wave wavelength .

[0032] Step 4.4: Calculate the global wave angular frequency. .

[0033] Calculate the discrete water depth values ​​for the adaptive window:

[0034] .

[0035] in, Indicates the first Discrete water depth values ​​at each water depth observation point, wave velocity , This is the acceleration due to gravity.

[0036] As described above, the main wavelength of the window Extraction is performed using the following steps:

[0037] The initial wave wavelength is calculated using the following formula. :

[0038] .

[0039] in, It is the acceleration due to gravity. For wave cycles, The initial water depth for GEBCO.

[0040] Set the initial size of the window to adapt. , This is an empirical coefficient.

[0041] For each beam, the position of the adaptive window slides continuously along the beam trajectory of the nearshore shallow water study area with a preset sliding step size.

[0042] For each adaptive window, perform the following iterative process independently at its position:

[0043] The reconstructed nearshore sea surface elevation sequence Perform Fast Fourier Transform (FFT) on the subsequence within the corresponding adaptive window to extract the wave wavelength corresponding to the wave with the largest amplitude. , Indicates the number of iterations.

[0044] like or Adjust the window size to fit the user. :

[0045] .

[0046] Return and perform the Fast Fourier Transform (FFT) again to extract the wave wavelength corresponding to the wave with the largest amplitude, until the extracted wave wavelength meets the physical constraints. The corresponding wave wavelength is denoted as the window dominant wavelength. The center of the adaptive window is the water depth observation point; This is the serial number of the water depth observation point.

[0047] As described above, step 4.3 includes the following steps:

[0048] The angle between the wave propagation direction and the ICESat-2 beam trajectory direction to be inverted is calculated based on the following formula:

[0049] .

[0050] in, The dominant wavenumber is the angle between the wave propagation direction and the beam trajectory direction. , The main wavelength; The phase difference between strong and weak beam signals; The spanning distance between strong and weak beams.

[0051] like True wave number Corrected wave wavelength ;like True wave number .

[0052] For nearshore shallow water research areas:

[0053] The dominant wavelength corresponding to each adaptive window The corresponding corrected wave wavelength is denoted as the corrected deep-water wave wavelength. .

[0054] For the entire deep water area:

[0055] Main wavelength The corresponding corrected wave wavelength is denoted as the corrected window dominant wavelength. .

[0056] As described above, step 5 specifically includes the following steps:

[0057] Step 5.1: Based on the selected nearshore shallow water study area, create a planar grid covering the nearshore shallow water study area; obtain the geographic coordinates of the coastline of the nearshore shallow water study area, and based on the preset offshore segment length... All water depth observation points and target grid points in the planar grid are divided into corresponding offshore distance segments.

[0058] Step 5.2: For each target grid point, using the target grid point as the center and the preset attenuation length parameter... Search for all water depth observation points within a circle of radius , and use them as neighborhood water depth observation points of the target grid point.

[0059] Step 5.3: For each target grid point, calculate its interpolation weights for all neighboring water depth observation points across all beams. The interpolation weights include spatial weights and segmentation weights.

[0060] For each target grid point, the water depth interpolation calculation is performed to obtain the corresponding interpolated water depth. The geographic coordinates of all target grid points and the corresponding interpolated water depths constitute the terrain digital elevation model (DEM). The water depth interpolation calculation for the target grid points specifically includes the following steps:

[0061] Calculate the interpolation weights of neighboring water depth observation points and interpolation weights Normalization is performed to obtain the normalized interpolation weights. ; For the first The first beam Spatial weights of neighborhood water depth observation points; For the first The first beam The segmented weights of each neighborhood depth observation point.

[0062] Calculate the interpolated water depth at the target grid points :

[0063] .

[0064] Indicates the first The first beam The water depth values ​​at each water depth observation point.

[0065] The spatial weights, as described above, are calculated based on the following formula:

[0066] .

[0067] in, For the first relative target grid point The first beam Spatial weights of neighborhood water depth observation points; and These represent the coordinate differences between the water depth observation point and the target grid point in the east-west and north-south directions of the plane coordinate system, respectively. This is the attenuation length parameter.

[0068] The segment weights are calculated based on the following formula:

[0069] .

[0070] in, For the first relative target grid point The first beam Segmented weights for each neighborhood depth observation point; and These are the numbers of the offshore distance segments where the neighboring water depth observation points are located and the numbers of the offshore distance segments where the target grid point is located; This is the attenuation factor.

[0071] A computer device includes a memory and a processor, the memory storing a computer program, the processor executing the computer program to implement steps 1 to 5 of the multi-temporal ICESat-2 turbid shallow water depth inversion method as described in any one of the claims.

[0072] Compared with the prior art, the present invention has the following beneficial effects:

[0073] (1) Robustness: The modeling of undulating sea surface is robust through gap sensing and dynamic buffering strategies. It is applicable to highly turbid waters, overcomes the limitations of optical methods, and can adapt to different sea conditions through multi-time series data fusion.

[0074] (2) High accuracy: On the one hand, the adaptive window fast Fourier transform (FFT) can effectively capture local changes in the wave spectrum and improve the accuracy of water depth inversion. On the other hand, the SDW fusion algorithm further optimizes the accuracy in the process of generating a continuous underwater digital elevation model (DEM) by introducing the core element of terrain change, "offshore distance", as a constraint.

[0075] (3) Continuity: The SDW fusion algorithm is used to achieve effective fusion of multi-temporal data and generate a complete and continuous underwater topographic digital elevation model (DEM), providing reliable basic data for marine engineering and ecological research.

[0076] The essential difference between this application and the prior art in the adaptive window fast Fourier transform (FFT) wavelength extraction process in step 4 is as follows:

[0077] Existing wave theory-based water depth inversion methods (whether using optical, SAR, or ICESat-2 data) generally employ Fast Fourier Transform (FFT) analysis with fixed or asymptotically long windows. This can contain mixed or incomplete wave signals, leading to inaccurate extraction of dominant swell features. In nearshore shallow waters, water depth varies dramatically, causing wave wavelengths to change significantly over short distances. Fixed windows cannot adapt to this variation: excessively long windows smooth out local details, introducing additional wavelength information that contaminates the shallow water signal; excessively short windows fail to capture complete waveform periods, rendering spectral analysis ineffective. This application proposes a physically constrained adaptive Fourier window strategy, directly addressing the inherent limitations of fixed windows and representing the key technical reason for achieving high-precision inversion. It ensures that the extracted wavelengths represent the local wave state, improving the accuracy of water depth inversion.

[0078] The essential difference between the offshore distance segmented dual-weight (SDW) fusion algorithm in step 5 of this application and existing technologies is as follows:

[0079] Existing methods invert data from a single satellite transit trajectory to obtain a discrete water depth profile distributed along the orbit. This result is linear, with its spatial coverage entirely dependent on the satellite's current flight trajectory, making it extremely sparse. While traditional spatial interpolation methods (such as inverse distance weighting and Kriging) can fuse results from multiple trajectories, these methods rely solely on the spatial geometry of the observation points. When fusing sparse, discrete water depth points like those from ICESat-2, they neglect a crucial prior knowledge about nearshore waters: water depth changes are highly correlated with offshore distance. Traditional methods, considering only spatial geometry, are prone to producing interpolation results that do not conform to actual terrain continuity in complex terrain areas near the coastline. This invention's pioneering Offshore Distance Segmented Dual Weight (SDW) fusion algorithm aims to fundamentally solve the aforementioned fusion challenge of "from sparse discrete points to a continuous and reasonable terrain surface." This scheme transcends the traditional interpolation paradigm based solely on spatial proximity, creatively introducing "offshore distance" as a core constraint variable controlling the trend of terrain change. The proposed solution is key to generating continuous and reasonable underwater topographic digital elevation models (DEMs). It solves the technical challenge of generating continuous surfaces from sparse points, making the fused results not only spatially smooth but also consistent with the nearshore topographic variation patterns, thereby significantly improving the accuracy of continuous products.

[0080] Furthermore, in step 3 of this application, the denoised photon data is reconstructed in profile:

[0081] ICESat-2 photon data often exhibits gaps of tens or even hundreds of meters along the track direction due to cloud cover, momentary instrument malfunctions, or failure of specular reflection from the sea surface. Traditional photon point cloud processing or profile reconstruction methods typically perform global curve fitting or filtering on the entire dataset. Treating this physically gapped data as a continuous whole for modeling would force the algorithm to fit a non-existent "sea surface," resulting in severely distorted and warped virtual profiles across the gaps. This introduces significant errors into subsequent wave spectrum analysis (FFT), fundamentally compromising the accuracy of water depth inversion.

[0082] This invention proposes a "gap-aware local B-spline fitting" method. This step outputs a high-quality, high-fidelity, continuous sea surface elevation sequence strictly based on the distribution of measured data for the entire inversion process. It preserves the undulation signal of the real waves to the maximum extent while eliminating systematic interference caused by missing data. Attached Figure Description

[0083] Figure 1 This diagram illustrates the overall process of the invention, showcasing the entire process of data preprocessing, adaptive window depth inversion, and multi-temporal fusion. The diagram clearly explains the logical relationships between the modules and the data flow. Here, FFT represents Fast Fourier Transform; GEBCO represents a world ocean seafloor topography map, denoted as GEBCO depth data.

[0084] Figure 2 This is a schematic diagram of sea surface reconstruction from a single ICESat-2 data point. (a) shows the result of conventional B-spline fitting, and (b) shows the result of gap-sensing local B-spline fitting. The horizontal axis represents the distance along the track, and the vertical axis represents the elevation. Yellow fitted points represent the reconstructed profile, gray points represent sea surface photons before denoising, and blue points represent sea surface photons after denoising. The unit is m (meter).

[0085] Figure 3 The diagram illustrates the inversion results of typical single-time-phase data from three ports, showcasing a comparison of the water depth inversion effects along a single ICESat-2 trajectory. The water depth inversion effects along a single ICESat-2 trajectory for (a) Port A, (b) Port C, and (c) Port B all include the following elements: blue scatter points represent sea surface photon points obtained after preprocessing and fitting; orange curves represent water depth results inverted using the adaptive Fourier window method of this invention; purple curves represent water depth results inverted using the traditional fixed Fourier window method; and green curves represent measured water depth data. RMSE is the root mean square error; R² is the coefficient of determination; and MAPE is the mean absolute percentage error.

[0086] Figure 4 This chart compares the accuracy of the adaptive Fourier window method and the fixed Fourier window method in the multi-time inversion results for three ports. In each sub-plot, the horizontal axis (x) represents the inverted water depth value, and the vertical axis (y) represents the measured water depth value. The red dashed line is the ideal fitting line, representing the theoretical reference that the inverted water depth is completely consistent with the measured water depth. The more concentrated the scatter points are near the red dashed line, the higher the inversion accuracy. Statistical indicators (RMSE: root mean square error; R²: coefficient of determination; MAPE: mean absolute percentage error) are attached to quantitatively compare the overall inversion accuracy of the two methods for ports A, C, and B. Specifically, (a) is a scatter plot of accuracy using the fixed window method for port A; (b) is a scatter plot of accuracy using the adaptive window method for port A; (c) is a scatter plot of accuracy using the fixed window method for port B; (d) is a scatter plot of accuracy using the adaptive window method for port B; (e) is a scatter plot of accuracy using the fixed window method for port C; and (f) is a scatter plot of accuracy using the adaptive window method for port C. (Unit: m) -1 Indicates rice -1 .

[0087] Figure 5 Scatter plots are provided to evaluate the accuracy of the fused continuous underwater topographic digital elevation model (DEM). In each subplot, the horizontal axis (x) represents the measured water depth, and the vertical axis (y) represents the inverted water depth. The red dashed line is the ideal fitting line, representing a theoretical reference where the inverted water depth is completely consistent with the measured water depth. The more concentrated the scatter points are near the red dashed line, the higher the inversion accuracy. Statistical indicators (Root Mean Square Error (RMSE), Mean Absolute Percentage Error (MAPE), and Coefficient of Decision (R²) are also provided. Specifically, (a) is a scatter plot evaluating the accuracy of the fused continuous underwater topographic digital elevation model (DEM) generated for the waters surrounding Port A, and (b) is a scatter plot evaluating the accuracy of the fused continuous underwater topographic digital elevation model (DEM) generated for the waters surrounding Port B. Detailed Implementation

[0088] To facilitate understanding and implementation of the present invention by those skilled in the art, the present invention will be further described in detail below with reference to embodiments. It should be understood that the embodiments described herein are for illustration and explanation only and are not intended to limit the present invention.

[0089] Example 1:

[0090] Multi-temporal ICESat-2 method for depth inversion in turbid shallow waters, such as Figure 1 As shown, the specific steps include:

[0091] Step 1: Identification and extraction of photon clusters on the sea surface.

[0092] A nearshore shallow water study area was selected, and multi-temporal ICESat-2 ATL03 level sea surface wave data covering the nearshore shallow water study area and the corresponding deep water area was acquired (in this embodiment, the area less than or equal to 40 kilometers from the coast is defined as the nearshore shallow water study area, and the area greater than 40 kilometers from the coast is defined as the deep water area). Raw photon point cloud data was extracted from the multi-temporal ICESat-2 ATL03 level sea surface wave data. The raw photon point cloud data includes photon point data corresponding to each photon under different beams, and the photon point data includes information such as the latitude, longitude, and elevation of the photons. The photons in the raw photon point cloud data were segmented along the orbital direction and partitioned vertically. For each segment along the orbital direction: based on the distribution of photon counts in the vertical partition, corresponding sea surface photon clusters were selected. The specific process of selecting sea surface photon clusters includes:

[0093] Step 1.1: Perform photon elevation truncation on the original photon point cloud data, retaining photons with elevation values ​​within a preset elevation range as initial selected photons. In this embodiment, photons with elevation values ​​within the range of [-50 meters, 50 meters] are retained. This range covers the possible range of photons on the sea surface, while excluding obvious outliers.

[0094] Step 1.2: For each beam in the original photon point cloud data, based on the beam trajectory, the photons along the track direction in the initial selection photons are segmented according to the preset track segment length (in this embodiment, the preset segment length is 2 kilometers), and each segment is recorded as a track direction segment.

[0095] Step 1.3: For each segment along the track direction: According to the preset vertical compartment height, the photons within the segment along the track direction are divided into vertical compartments, each of which is denoted as a vertical compartment. In this embodiment, a vertical compartment with a height of 1 meter is used, and a photon elevation histogram is generated to represent the number of photons corresponding to each vertical compartment. This step initially selects the photon elevation range of the sea surface photon cluster by analyzing the distribution of photons within the vertical compartments.

[0096] Step 1.4: To adapt to wave fluctuations and noise effects, a dynamic buffer strategy is adopted. For each segment along the track direction, based on the photon count of the vertical bin with the highest photon count, and the photon count in the vertical bin adjacent to the vertical bin with the highest photon count, an upper buffer and a lower buffer are constructed according to the following formula. The elevation extraction range of the corresponding sea surface photon clusters is adaptively set. .

[0097] in, This represents the height of the upper boundary of the vertical bin with the highest photon count. The lower boundary height of the vertical bin with the highest photon count. This indicates the height difference of the lower buffer zone. Indicates the height difference of the upper buffer zone:

[0098] .

[0099] .

[0100] in, The vertical compartment height is the height specified in this embodiment. The value is 1 meter; This represents the number of photons in the vertical bin with the largest number of photons. This represents the number of photons in the first vertical bin immediately below the vertical bin with the highest number of photons (at a lower photon elevation). This represents the number of photons in the first vertical bin immediately above the vertical bin with the highest number of photons (higher photon elevation).

[0101] Number of photons and This reflects the intensity of photons diffusing upwards and downwards from the sea surface. The ratio ( )and( This determines the number of "expansion steps" the buffer extends downwards and upwards (each expansion step equals the vertical bin height). ).ratio( )and( The larger the value, the stronger the signal in the vertical bin adjacent to the vertical bin with the largest number of photons. The buffer then expands further in that direction, thus adaptively enveloping the photons on the undulating sea surface caused by waves.

[0102] The upper boundary of the vertical bin with the largest number of photons and lower boundary Based on this, the elevation values ​​of each segment along the track fall within the corresponding interval. Photons within the ocean are considered as surface photons, and the collection of surface photons is called a surface photon cluster.

[0103] This buffer zone design can effectively adapt to wave fluctuations under different sea conditions, ensuring the complete extraction of photons from the sea surface.

[0104] Step 2: Denoise the photon clusters on the sea surface.

[0105] The sea surface photon clusters extracted in Step 1 still contain a large number of noisy photons, requiring further denoising. This invention employs a two-step denoising method: sequentially using a density-based spatial clustering (DBSCAN) denoising algorithm and a two-dimensional sliding window mid-range filter for two-step denoising. The specific process is as follows:

[0106] Step 2.1: The first step of denoising is performed on the sea surface photons in the sea surface photon cluster using the density-based spatial clustering (DBSCAN) algorithm. The specific parameter set in this embodiment is: denoising neighborhood radius. Meters, minimum number of samples for noise reduction .

[0107] Step 2.2: A two-dimensional sliding window mid-range filter is used for the second step of fine denoising. The length of the denoising window and the denoising sliding step size are taken as the diameter of the ICESat-2 spot (approximately 17 meters). For each denoising window sliding horizontally, the center point in the vertical direction (height direction) of the denoising window is calculated based on the median elevation of all sea surface photons covered by the length of the denoising window. Then, using the height of the center point in the vertical direction of the denoising window as a reference, the window is extended upwards and downwards by a preset half-height (in this embodiment, the preset half-height is 0.5 meters), thereby forming a vertical screening range with a window height of 1 meter. In this embodiment, the sliding step size is set to 17 meters, and the center point of the denoising window is traversed to iterate through the sea surface photon clusters. By reading the sc_orient (satellite ascent / descent orbit) parameter from the multi-temporal ICESat-2 ATL03 level sea surface wave data file, beams are distinguished into strong and weak beams. When the sc_orient value is 1, beams gt1r, gt2r, and gt3r are identified as strong beams, while beams gt1l, gt2l, and gt3l are identified as weak beams. When the sc_orient value is 0, beams gt1r, gt2r, and gt3r are identified as weak beams, while beams gt1l, gt2l, and gt3l are identified as strong beams. The minimum sample number parameter of the two-dimensional sliding window value filter is dynamically adjusted based on the beam strength, and the minimum sample number in the strong beam neighborhood is determined. Minimum number of samples for weak beam .

[0108] This step further filters out residual local noise, ultimately outputting high-quality signal photons that can be used for sea surface profile reconstruction.

[0109] Step 3: Reconstruct the wave surface profile.

[0110] The distribution of the denoised sea surface photon cluster along the orbital direction is still irregular and discrete, requiring the reconstruction of a segmented, continuous sea surface profile. This invention utilizes a gap-aware local B-spline fitting algorithm to reconstruct a continuous, smooth, undulating sea surface profile, generating a high-quality sea surface high-order sequence. In this embodiment, the space where the orbital distance between adjacent sea surface photons in the denoised photon cluster is greater than a preset gap threshold is used as the gap. Based on the gap, the sea surface photons corresponding to each beam are divided into multiple continuous segments. For each continuous segment, based on the orbital distance and elevation of the corresponding sea surface photon, independent smooth B-spline fitting and resampling are performed to generate the reconstructed sea surface high-order sequence. The specific process is as follows:

[0111] Step 3.1: For each beam, based on the beam trajectory, calculate the along-track spacing between adjacent sea surface photons in the denoised sea surface photon cluster. Spaces where the along-track spacing between adjacent sea surface photons is greater than a preset gap threshold (50 meters in this embodiment) are considered gaps. These gaps are usually caused by cloud cover or instrument malfunction and require special handling; see steps 3.2-3.3 for details.

[0112] Step 3.2: Using the gaps identified in Step 3.1 as boundaries, divide the sea surface photons in the denoised sea surface photon cluster from Step 2 into multiple continuous segments; for each continuous segment, based on the orbital distance and elevation of the corresponding sea surface photon, independently perform smooth B-spline fitting, generating a smooth and continuous B-spline fitting curve function, such as... Figure 2 As shown, this curve describes the variation of sea surface photon elevation with distance along the track within a corresponding continuous segment.

[0113] The "along-track distance" is a continuously accumulated physical quantity within the photon point cloud data of each beam, used to precisely describe the absolute position of a photon on its respective beam trajectory. Its calculation does not rely on arbitrarily chosen "reference points," but is based on the segment geometry information provided in the ICESat-2 ATL03 data product. For each photon, its along-track distance... ,in The reference distance along the orbit from the starting point of the photon's orbital segment. This represents the offset of the photon along the orbit relative to the starting point of the corresponding orbital segment. This refers to the inherent "track segment" number provided by the product, used to calculate the absolute track distance of the photon (not the "track direction segment" set in step 1 of this invention for signal extraction).

[0114] This method calculates a set of global, segmented, continuous, and physically meaningful along-track distance coordinates for all sea surface photons within the same beam, serving as the spatial basis for profile reconstruction. In subsequent processing, even if the photon data is divided into multiple continuous segments due to gaps, photons within each continuous segment still use this unified set of along-track distance coordinates, thus ensuring that the profiles fitted to each continuous segment of the same beam are accurately aligned in spatial position. In this embodiment, during smooth B-spline fitting, the smoothing factor corresponding to the strong beam... Smoothing factor corresponding to weak beam ,in For the first The number of photons in a continuous segment.

[0115] Step 3.3: For each beam: For each continuous segment corresponding to the beam, based on the number of evaluation points... The corresponding B-spline fitting curves are resampled to obtain the corresponding reconstructed sub-sea surface height sequences; all reconstructed sub-sea surface height sequences are then spliced ​​together to obtain the reconstructed sea surface height sequence. .

[0116] After obtaining the B-spline fitting curve (a mathematical function) for each continuous segment, in order to obtain a set of regularly sampled data sequences that facilitate subsequent analysis (Fourier Transform, FFT), it is necessary to sample the data within the domain (distance range along the track) of the B-spline fitting curve at a certain resolution (i.e., the number of evaluation points). Resampling is performed to calculate a new series of regularly spaced points (distance along the track, elevation). This sequence, composed of resampled points from all consecutive segments, is the reconstructed sea surface elevation sequence. :

[0117] .

[0118] The reconstructed sea surface height sequence, which represents the segmented continuous sea surface height sequence reconstructed from the fitted data along the entire beam trajectory, is the basic data for subsequent wave spectrum analysis. Index representing a continuous segment , Indicates the number of consecutive segments on a beam trajectory; Representing the The orbital distance and geographic coordinates of the resampled profile points within a continuous segment. To fit the sea surface elevation value, characterize and The corresponding sea surface elevation value is obtained from the B-spline fitting curve. It represents the union of sets.

[0119] The core innovation of step 3 lies in the strategies of "gap sensing" and "segmented independent fitting," as well as the detail of adaptively setting the smoothing factor based on the beam type and the number of photons within a segment. These designs collectively solve the problem of severe distortion and spurious signals generated by traditional fitting methods when there are large gaps in the data.

[0120] Multi-temporal ICESat-2 ATL03-level sea surface wave data from the nearshore shallow water research area and adjacent deep water area were processed using steps 1 to 3 (sea surface photon extraction, denoising, and wave profile reconstruction) to ensure data quality and comparability.

[0121] Step 4: Invert water depth based on linear dispersion relation.

[0122] Based on the reconstructed sea surface elevation sequence Using Fast Fourier Transform (FFT), the wave wavelengths corresponding to the waves with the largest amplitudes were extracted and corrected in both the nearshore shallow water and deep water study areas. Based on the corrected wave wavelengths, discrete water depth values ​​at different observation points in the nearshore shallow water study area were then extracted. The specific process is as follows:

[0123] Step 4.1: Follow the reconstructed sea surface elevation sequence The corresponding resampling profile points are located in nearshore shallow water or deep water areas, which will reconstruct the sea surface elevation sequence. Divided into reconstructed nearshore sea surface elevation series and the reconstruction of deep-water sea surface elevation sequence .

[0124] Step 4.2: Based on Fast Fourier Transform (FFT), extract the wave wavelengths corresponding to the waves with the largest amplitudes in the nearshore shallow water study area and the deep water area, respectively.

[0125] (1) For deep water areas:

[0126] Corresponding to the entire deep water area, the corresponding reconstructed deep water surface elevation sequence for each beam. Perform a Fast Fourier Transform (FFT) to extract the wave wavelength corresponding to the wave with the maximum amplitude, and denote it as the main wave wavelength in the deep water area. .

[0127] (2) For nearshore shallow water study areas:

[0128] Based on the prior initial water depth, an adaptive window is initialized; for each beam: the adaptive window slides along the beam trajectory of the nearshore shallow water study area based on a preset sliding step size; the reconstructed nearshore sea surface elevation is then programmed. The size of the adaptive window is iteratively adjusted using Fast Fourier Transform (FFT) and a preset wave wavelength range for the subsequences corresponding to each location of the adaptive window, and the main wave wavelength within the adjusted adaptive window range is extracted. This allows for the adaptive extraction of the local dominant wavelength. The center of the adaptive window corresponds to the water depth observation point. This is the sequence number of the water depth observation point; each water depth observation point corresponds to an adaptive window at a specific location.

[0129] First, initialize the adaptive window: use the existing publicly available GEBCO water depth data as the initial water depth, and calculate the initial wave wavelength using linear wave theory. :

[0130] .

[0131] in, It is the acceleration due to gravity. For wave cycles, The initial water depth for GEBCO is given; in this embodiment, the Brent numerical iteration method is used to solve the implicit equations above.

[0132] Set the initial size of the window to adapt. ,in This is an empirical coefficient. .

[0133] Then, for each beam, the position of the adaptive window continuously slides along the beam trajectory of the nearshore shallow water study area with a preset sliding step size. In this embodiment, the sliding step size is... .

[0134] For each adaptive window, perform the following iterative process independently at its position:

[0135] The reconstructed nearshore sea surface elevation sequence Perform Fast Fourier Transform (FFT) on the subsequence within the corresponding adaptive window to extract the wave wavelength corresponding to the wave with the largest amplitude. , Indicates the number of iterations.

[0136] like or Adjust the window size to fit the user. (Unit: meters):

[0137] .

[0138] Return to the previous step and perform a Fast Fourier Transform (FFT) again to extract the wave wavelength corresponding to the wave with the largest amplitude. Repeat this process until the extracted wave wavelength meets the physical constraints. The corresponding wave wavelength is denoted as the window dominant wavelength. The center of the adaptive window is the water depth observation point; This is the sequence number of the water depth observation point; each water depth observation point corresponds to an adaptive window at a specific location.

[0139] Step 4.3: Wavelength correction and wave direction correction.

[0140] By utilizing the observation phase difference between strong and weak beams, the angle between the wave propagation direction and the beam trajectory direction is corrected, thereby correcting the corresponding dominant wave wavelength. The dominant wave wavelength includes the dominant wave wavelength in deep water. and the main wavelength of the window The corrected deep-water wave wavelengths were obtained accordingly. and the corrected window main wavelength .

[0141] The angle between the wave propagation direction and the ICESat-2 beam trajectory direction to be inverted is calculated based on the following formula:

[0142] .

[0143] in, The dominant wavenumber is the angle between the wave propagation direction and the beam trajectory direction. , The main wavelength; The phase difference between strong and weak beam signals; The inter-track spacing between the strong and weak beams is known to be 90m; the phase difference... and track spacing All are sea surface elevation sequence reconstructed from step 3 get.

[0144] like Based on the included angle Correction yields the true wave number. Thus, the corrected wave wavelength is obtained. ;like Negligible, actual wave number .

[0145] For nearshore shallow water research areas:

[0146] The dominant wavelength corresponding to each adaptive window The corresponding corrected wave wavelength is denoted as the corrected window dominant wavelength. .

[0147] For the entire deep water area:

[0148] Main wavelength The corresponding corrected wave wavelength is denoted as the corrected deep-water wave wavelength. .

[0149] Step 4.4: Based on the deep-water dispersion relation in linear wave theory, calculate the global wave angular frequency representing the incident wave state. .

[0150] The global wave angular frequency obtained in the deep water area Corrected window main wavelength corresponding to the nearshore shallow water study area By combining these methods and substituting the linear dispersion relation, the adaptive windows for the nearshore shallow water study area are obtained through inversion calculation. Discrete water depth values ​​within:

[0151] .

[0152] in, Indicates the first Discrete water depth values ​​at each water depth observation point The serial number of the water depth observation point; wave velocity .

[0153] This invention independently executes the entire processing flow from step 1.1 to step 4.4 for each beam. Each beam generates a set of independent discrete water depth values ​​for water depth observation points located at different adaptive window centers within its own nearshore shallow water study area. The water depth calculation described in step 4.4 is performed separately for each effective adaptive window of each beam. The discrete water depth value of a "water depth observation point" specifically refers to the inverted water depth value of a certain beam at a certain adaptive window position. These discrete water depth values ​​generated by all beams will be used together as the input dataset for step 5 (multi-temporal data fusion).

[0154] Step 5: Generate a spatially continuous underwater terrain digital elevation model (DEM).

[0155] A planar grid covering the nearshore shallow water study area was created. The offshore distances of all depth observation points and target grid points within the planar grid were calculated and segmented. Within a circle centered on each target grid point and with a decay length parameter as the radius, the spatial weights and segmentation weights of neighboring depth observation points were comprehensively calculated to obtain the interpolation weights for these neighboring points. Based on these interpolation weights and discrete depth values, depth interpolation was performed on all target grid points within the nearshore shallow water study area to generate an underwater topographic digital elevation model (DEM). The specific process is as follows:

[0156] Step 5.1: Based on the selected nearshore shallow water study area, create a planar grid covering the nearshore shallow water study area; obtain the geographic coordinates of the coastline of the nearshore shallow water study area, and based on the preset offshore segment length... All water depth observation points and target grid points in the planar grid are divided into corresponding offshore distance segments.

[0157] Step 5.1.1: Create a planar grid covering the nearshore shallow water study area, and obtain the geographic coordinates of the coastline of the nearshore shallow water study area:

[0158] Based on the nearshore shallow water study area selected in step 1, a regular planar grid covering the nearshore shallow water study area is created according to the specified spatial resolution. Each grid intersection is a target grid point to be interpolated. This embodiment uses a spatial resolution of 30 meters and utilizes a TIFF image to create a planar grid covering the nearshore shallow water study area. Each target grid point has a unique corresponding geographic coordinate.

[0159] Coastline vectorization: Based on the World Imagery Wayback remote sensing image base map provided by ArcGIS, coastline vector data is generated, including the geographic coordinates of each point on the coastline. Its core purpose is to calculate the accurate offshore distance for all target grid points and water depth observation points during subsequent multi-temporal data fusion. This is a key input for implementing the offshore distance segmentation and offshore distance segmentation dual-weight fusion algorithm (denoted as the SDW algorithm). The nearshore shallow water study area in step 1 is delineated, covering only a few kilometers, far from the deep water area. In step 4 (single-temporal water depth inversion), this high-precision coastline was not relied upon when distinguishing between the nearshore and deep water areas.

[0160] Step 5.1.2: Calculate and segment the offshore distances for multi-temporal discrete water depth points, so that all water depth observation points and target grid points are assigned to the corresponding offshore distance segments:

[0161] Calculate the offshore distance from all water depth observation points and target grid points to the coastline, and divide the data into segments: for each discrete water depth value obtained from the adaptive window water depth inversion in step 4... Calculate the shortest Euclidean distance between the corresponding water depth observation point and the coastline. The shortest Euclidean distance reflects the distance of the water depth observation point from the shore and is the basis for segmentation. Calculate the segment length based on the preset distance from the shore. The offshore distance range is divided into multiple offshore distance calculation segments. This range effectively balances the precision required to capture terrain changes with the number of samples within each segment.

[0162] Step 5.2: For each target grid point, taking the target grid point as the center, within the distance attenuation influence range determined by the spatial weight attenuation length (i.e., using the preset attenuation length parameter)... Search all depth observation points within a circle of radius (where the circle is the radius of the target grid point) and use them as neighborhood depth observation points.

[0163] Step 5.3: Employ the Offshore Distance Segmented Dual-Weight Fusion Algorithm (SDW algorithm) to calculate the spatial weight and segmented weight, and then fuse the spatial weight and segmented weight to calculate the interpolated water depth at the target grid points. The specific process is as follows:

[0164] Step 5.3.1: For each target grid point of the rule that needs to generate an underwater terrain digital elevation model (DEM), calculate its interpolation weights for the neighboring water depth observation points in all beams. The interpolation weights include spatial weights and segment weights.

[0165] (1) Spatial weights:

[0166] Spatial weights are calculated using a Gaussian decay function to measure the planar spatial proximity between neighboring depth observation points and the target grid point. The weighting relative to the target grid point is calculated based on the following formula. The first beam Spatial weight of each neighborhood depth observation point :

[0167] .

[0168] in, and These represent the coordinate differences in the east-west and north-south directions between the water depth observation point (i.e., the location of the discrete water depth values ​​obtained in step 4) and the target grid point in a plane coordinate system (such as UTM coordinates). The attenuation length parameter (used to control the spatial weight attenuation rate, preset in this embodiment) ).

[0169] (2) Segment weights:

[0170] The segmented weight calculation is based on the difference in offshore distance between the target grid point and its neighboring water depth observation points. The weight is calculated relative to the target grid point using the following formula. The first beam Piecewise weights of neighborhood depth observation points :

[0171] .

[0172] in, and These are the numbers of the offshore distance segments where the neighboring water depth observation points are located and the numbers of the offshore distance segments where the target grid point is located; This is the decay factor (usually 0.5), which controls the decay rate of the segment weights.

[0173] Step 5.3.2: For each target grid point, perform water depth interpolation calculation to obtain the corresponding interpolated water depth. The geographic coordinates of all target grid points and their corresponding interpolated water depths constitute a spatially continuous underwater topographic digital elevation model (DEM). The water depth interpolation calculation for the target grid points specifically includes the following steps:

[0174] Calculate the interpolation weights of neighboring water depth observation points and interpolation weights Normalization is performed to obtain the normalized interpolation weights. .

[0175] The interpolated water depth at the target grid points is calculated using a weighted average. :

[0176] .

[0177] Interpolation depth of target grid points This represents the final calculated water depth value at a specific grid location on the underwater terrain digital elevation model (DEM) to be generated. Indicates the first The first beam The water depth values ​​at each water depth observation point are obtained from step 4.

[0178] Example 2:

[0179] To demonstrate the effectiveness of this invention, the multi-temporal ICESat-2 method for inverting the depth of turbid shallow waters was used. Three typical turbid shallow water areas—Port A, Port B, and Port C—were selected for inversion tests using an adaptive Fourier window. The multi-temporal results from the inversions of Port A and Port B were then fused to generate a continuous underwater topographic digital elevation model (DEM). Table 1 shows the parameter settings.

[0180] Table 1 Parameter Setting Instructions

[0181]

[0182] This invention provides more accurate water depth inversion results. By taking a single typical trajectory data point, the inversion accuracy of the fixed window method and the adaptive window method of this invention are compared. The inversion results of typical single-time-phase data from three ports are as follows: Figure 3As shown in Table 2, typical results for the three regions are compared, demonstrating a significant improvement in accuracy using the adaptive window method. Here, RMSE represents the root mean square error; R² is the coefficient of determination; and MAPE is the mean absolute percentage error.

[0183] Table 2. Accuracy Evaluation of Inversion Results Based on Typical Single-Phase Data of the Invention

[0184]

[0185] All available ICESat-2 transit data for the study area from October 2018 to October 2024 were collected, and inversion was performed using two different methods. A comparison of the accuracy of the adaptive Fourier window and fixed Fourier window methods in the multi-temporal inversion results for the three ports is presented below. Figure 4 As shown in Table 3, the overall accuracy comparison of the three typical regions is presented. The method of this invention exhibits stable and higher accuracy in different regions.

[0186] Table 3. Accuracy Evaluation of Inversion Results Based on Multi-Temporal Data in this Invention

[0187]

[0188] The ICESat-2 trajectories of port areas A and B are consistent with the swell direction, resulting in abundant effective wave observation data. For all adaptive window inversion depth points in these two areas, the SDW algorithm proposed in this invention is used to fuse them to generate a 30-meter resolution underwater topographic digital elevation model (DEM). This model is then compared with measured water surface data, showing varying degrees of improvement compared to the original inversion results. The multi-time inversion results from the two ports are fused using SDW to generate a continuous underwater topographic digital elevation model (DEM). The scatter plot evaluating the accuracy of the fused continuous underwater topographic digital elevation model (DEM) is shown below. Figure 5 As shown in Table 4, the accuracy assessment results show that the fused continuous terrain digital elevation model (DEM) matches the measured data well, with an overall accuracy of up to 2 meters.

[0189] Table 4. Accuracy Evaluation Table of the Invention for Multi-Temporal Results Fusion of Continuous Terrain

[0190]

[0191] Example 3:

[0192] The multi-temporal ICESat-2 turbid shallow water depth inversion apparatus is used to implement the multi-temporal ICESat-2 turbid shallow water depth inversion method described in Example 1, including:

[0193] The sea surface photon cluster extraction module is used to implement step 1 of embodiment 1.

[0194] The noise reduction module is used to implement step 2 of embodiment 1.

[0195] The sea surface elevation sequence reconstruction module is used to implement step 3 in embodiment 1.

[0196] The discrete water depth calculation module is used to implement step 4 in Example 1.

[0197] The DEM generation module is used to implement step 5 in Example 1.

[0198] Example 4:

[0199] In this embodiment, a computer device is also provided, including a memory and a processor. The memory stores a computer program, and the processor executes the computer program to implement the steps of the multi-temporal ICESat-2 turbid shallow water depth inversion method described in Embodiment 1 above.

[0200] Example 5:

[0201] In this embodiment, a computer-readable storage medium is provided, on which a computer program is stored. When the computer program is executed by a processor, it implements the steps of the multi-temporal ICESat-2 turbid shallow water depth inversion method described in Embodiment 1 above.

[0202] Example 6:

[0203] In this embodiment, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps of the multi-temporal ICESat-2 turbid shallow water depth inversion method described in Embodiment 1 above.

[0204] It should be noted that the specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains can make various modifications or additions to the described specific embodiments or use similar methods to substitute them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A multi-temporal ICESat-2 method for depth inversion in turbid shallow waters, characterized in that, Includes the following steps: Step 1: Acquire multi-temporal ICESat-2 ATL03 level sea surface wave data covering the nearshore shallow water research area and the corresponding deep water area; extract raw photon point cloud data from the multi-temporal ICESat-2 ATL03 level sea surface wave data, and segment the photons in the raw photon point cloud data by the track direction and the vertical direction; for each track direction segment: based on the distribution of the number of photons in the vertical direction segment, select the corresponding sea surface photon clusters; Step 2: Denoise the photon clusters on the sea surface; Step 3: The space between adjacent sea surface photons in the denoised sea surface photon cluster that is greater than a preset gap threshold is used as the gap. Based on the gap, the sea surface photons corresponding to each beam are divided into multiple continuous segments. For each continuous segment, based on the along-track distance and elevation of the corresponding sea surface photon, smooth B-spline fitting and resampling are performed independently to generate a reconstructed sea surface elevation sequence. ; Step 4: Based on the reconstructed sea surface elevation sequence Using Fast Fourier Transform (FFT), the wave wavelengths corresponding to the waves with the largest amplitudes are extracted and corrected in the nearshore shallow water study area and the deep water area, respectively. Based on the corrected wave wavelengths, discrete water depth values ​​of different water depth observation points in the nearshore shallow water study area are extracted. Step 5: Create a planar grid covering the nearshore shallow water study area. Calculate and segment the offshore distances for all depth observation points and target grid points within the planar grid, dividing the offshore distance range into multiple offshore distance calculation segments. Calculate the interpolation weights for neighboring depth observation points centered on each target grid point. The interpolation weights include spatial weights and segment weights. Spatial weights measure the planar spatial proximity between neighboring depth observation points and target grid points. Segment weights calculate the difference between the offshore distance segment where the target grid point is located and the offshore distance segment where the neighboring depth observation points are located. Based on the interpolation weights of the neighboring depth observation points and the discrete depth values, perform depth interpolation calculations for all target grid points to generate an underwater topographic digital elevation model (DEM).

2. The method for depth inversion of turbid shallow waters using multi-temporal ICESat-2 as described in claim 1, characterized in that, The screening of the sea surface photon clusters specifically includes the following steps: Step 1.1: Extract photon elevations from the original photon point cloud data, retaining photons whose elevation values ​​are within a preset elevation range, and use them as initial photons; Step 1.2: For each beam in the original photon point cloud data, based on the beam trajectory, the photons along the track direction in the initial selection photons are segmented according to the preset track segment length, and each segment is recorded as a track direction segment. Step 1.3: For each segment along the track direction: According to the preset vertical box height, the photons in the segment along the track direction are divided into vertical boxes. Each vertical box is called a vertical box. Step 1.4: For each segment along the track direction, adaptively set the corresponding elevation extraction range of the sea surface photon cluster. ; This represents the height of the upper boundary of the vertical bin with the highest photon count. The lower boundary height of the vertical bin with the highest photon count. This indicates the height difference of the lower buffer zone. Indicates the height difference of the upper buffer zone: ， ， in, This refers to the vertical compartment height; This represents the number of photons in the vertical bin with the largest number of photons. This represents the number of photons in the first vertical bin immediately below the vertical bin with the largest number of photons. This represents the number of photons in the first vertical box immediately above the vertical box with the largest number of photons. The elevation values ​​of each segment along the track direction fall within the corresponding intervals. Photons within the ocean are considered as surface photons, and the collection of surface photons is called a surface photon cluster.

3. The method for inverting the depth of turbid shallow waters using multi-temporal ICESat-2 as described in claim 1, characterized in that, Step 2 uses a density-based spatial clustering denoising algorithm and a two-dimensional sliding window midpoint filter to denoise the sea surface photons in the sea surface photon cluster.

4. The method for inverting the depth of turbid shallow waters using multi-temporal ICESat-2 as described in claim 1, characterized in that, Step 3 specifically includes the following steps: Step 3.1: For each beam, based on the beam trajectory, calculate the along-track spacing between adjacent sea surface photons in the denoised sea surface photon cluster, and use the space where the along-track spacing is greater than the preset gap threshold as the gap. Step 3.2: Using the gaps identified in Step 3.1 as boundaries, divide the sea surface photons in the denoised sea surface photon cluster in Step 2 into multiple continuous segments; for each continuous segment, perform smooth B-spline fitting independently based on the orbital distance of the corresponding sea surface photon and the elevation of the sea surface photon to generate the corresponding B-spline fitting curve function. Step 3.3: For each beam: For each continuous segment corresponding to the beam, based on the number of evaluation points... , For the first The number of photons within a continuous segment is used to resample the corresponding B-spline fitting curve, resulting in the reconstructed sub-sea surface height sequence. All reconstructed sub-sea surface height sequences are then concatenated to obtain the reconstructed sea surface height sequence. .

5. The method for depth inversion of turbid shallow waters using multi-temporal ICESat-2 as described in claim 1, characterized in that, Step 4 specifically includes the following steps: Step 4.1: Follow the reconstructed sea surface elevation sequence The corresponding resampling profile points are located in nearshore shallow water or deep water areas, which will reconstruct the sea surface elevation sequence. Divided into reconstructed nearshore sea surface elevation series and the reconstruction of deep-water sea surface elevation sequence ; Step 4.2: Serialize the reconstructed deep-water sea surface elevation for each beam. Perform a Fast Fourier Transform (FFT) to extract the wave wavelength corresponding to the wave with the maximum amplitude, and denote it as the main wave wavelength in the deep water area. ; Based on the prior initial water depth, an adaptive window is initialized; for each beam: the adaptive window slides along the beam trajectory of the nearshore shallow water study area based on a preset sliding step size; the reconstructed nearshore sea surface elevation is then programmed. The size of the adaptive window is iteratively adjusted using Fast Fourier Transform (FFT) and a preset wave wavelength range for the subsequences corresponding to each location of the adaptive window, and the main wave wavelength within the adjusted adaptive window range is extracted. The center of the adaptive window is the water depth observation point; This refers to the serial number of the water depth observation point; Step 4.3: Use the observation phase difference between strong and weak beams to correct the angle between the wave propagation direction and the beam trajectory direction, thereby correcting the corresponding dominant wave wavelength. The dominant wave wavelength includes the dominant wave wavelength in deep water. and the main wavelength of the window The corrected deep-water wave wavelengths were obtained accordingly. and the corrected window main wavelength ; Step 4.4: Calculate the global wave angular frequency. ; Calculate the discrete water depth values ​​for the adaptive window: ， in, Indicates the first Discrete water depth values ​​at each water depth observation point, wave velocity , This is the acceleration due to gravity.

6. The method for depth inversion of turbid shallow waters using multi-temporal ICESat-2 according to claim 5, characterized in that, The main wavelength of the window Extraction is performed using the following steps: The initial wave wavelength is calculated using the following formula. : ， in, It is the acceleration due to gravity. For wave cycles, The initial water depth for GEBCO; Set the initial size of the window to adapt. , This is an empirical coefficient; For each beam, the position of the adaptive window slides continuously along the beam trajectory of the nearshore shallow water study area with a preset sliding step size; For each adaptive window, perform the following iterative process independently at its position: The reconstructed nearshore sea surface elevation sequence Perform Fast Fourier Transform (FFT) on the subsequence within the corresponding adaptive window to extract the wave wavelength corresponding to the wave with the largest amplitude. , Indicates the number of iterations; like or Adjust the window size to fit the user. : ， Return and perform the Fast Fourier Transform (FFT) again to extract the wave wavelength corresponding to the wave with the largest amplitude, until the extracted wave wavelength meets the physical constraints. The corresponding wave wavelength is denoted as the window dominant wavelength. The center of the adaptive window is the water depth observation point; This is the serial number of the water depth observation point.

7. The method for depth inversion of turbid shallow waters using multi-temporal ICESat-2 according to claim 5, characterized in that, Step 4.3 includes the following steps: The angle between the wave propagation direction and the ICESat-2 beam trajectory direction to be inverted is calculated based on the following formula: ， in, The dominant wavenumber is the angle between the wave propagation direction and the beam trajectory direction. , The main wavelength; The phase difference between strong and weak beam signals; The inter-track spacing between strong and weak beams; like True wave number Corrected wave wavelength ;like True wave number ; For nearshore shallow water research areas: The dominant wavelength corresponding to each adaptive window The corresponding corrected wave wavelength is denoted as the corrected window dominant wavelength. ; For the entire deep water area: Main wavelength The corresponding corrected wave wavelength is denoted as the corrected deep-water wave wavelength. .

8. A computer device comprising a memory and a processor, the memory storing a computer program, the processor executing the computer program to implement steps 1 to 5 of the multi-temporal ICESat-2 turbid shallow water depth inversion method according to any one of claims 1 to 7.

Images