A new extraction and grading method for underwater topographic photons of photon-counting LiDAR with different observation conditions

ABSTRACT

Spaceborne photon-counting light detection and ranging (LiDAR) have been extensively applied in shallow-water bathymetry. The density of underwater topographic photons (UTP) varies and is discontinuous due to sunlight noise, beam intensity, and seabed reflectivity, which differ from the land photon distribution due to the attenuation of water. Therefore, a general method for extracting and grading UTP is still lacking. We propose an active contour method combined with a variable convolution kernel method to calculate the photon range by considering the energy contributions of adjacent photons. Adaptive parameters under different observation conditions were determined to obtain the optimal convolution kernel using a kernel ridge regression model. This implies that the number of photons contained in the buffer zone was largest after the extracted UTP was fitted to a curve. Quantitative and qualitative verifications proved that the method performed well under different conditions. The photons obtained by the energy functional and the curve obtained by the fitting method were then used to grade the photons. Finally, an online developed UTP dataset and extraction framework were proposed to provide an applicable method for current and subsequent spaceborne photon-counting LiDAR.

KEYWORDS:

• Photon-counting LiDAR
• photons denoising
• underwater topographic photons
• active contours
• kernel ridge regression
• ICESat-2

1. Introduction

Information on shallow water depth is significant in the marine geodesy and cartography fields, mainly for shallow water in the coastal zone and shallow water near island reefs. Shallow water depth provides basic data for fishery culture, coastal zone management, and ship navigation (Bergsma et al. 2021). Coral reef ecosystems, the most unique ecosystems on Earth, are found in shallow waters. Therefore, studying the acquisition of shallow water depth is of great importance (Dong et al. 2019). Shallow water-sounding techniques for coastal zones and island reefs are diverse, but marine surveying and mapping still face challenges in obtaining large-range, high-resolution, and high-accuracy shallow water depth (IHO 2018). Various observation technologies have been used for shallow water bathymetry, but measured data on shallow water depth are still lacking, despite the development of satellite remote sensing technology (Li et al. 2023; Wölfl et al. 2019). Although ICESat-2 has been used for bathymetry, the UTP extraction of ATL03 still lacks a universal algorithm.

Shallow water depth acquisition can be divided into active and passive inversions. Active detection primarily includes acoustic and laser ranging (Costa, Battista, and Pittman 2009). Shipborne sonar (single-beam, multi-beam, etc.) is expensive and time-consuming, making it difficult to complete water depth mapping in shallow and sensitive areas near reefs and rocks (Su et al. 2018). Airborne lidar is an effective method for shallow sea sounding in coastal zones, which is not suitable for large-scale observations, and it is difficult to reach the ocean area (Irish and Lillycrop 1999). Passive water depth inversion includes empirical, physical, and semianalytical models. Passive inversion is far less accurate than active detection. However, its advantage is that it can be observed over a large area and long time series and its cost is low (Kutser et al. 2020). Research on active and passive fusion sounding is ongoing. In addition, SAR can estimate the depth of shallow water by using the change in the detected sea surface roughness. However, its mechanism is complex, and the requirements for wind speed and velocity are strict (Pereira et al. 2019). Dual-media photogrammetry can be used in limited areas with calm water surfaces, high water transparency, and rich bottom texture by exploiting the relative geometric relationship of optical stereo-image pairs (Hodúl et al. 2018). However, the application scenarios for both methods are very limited.

In recent years, massive amounts of image data have promoted the development of active and passive water depth inversions. Empirical and semi-analytical models depend heavily on the water depth label (Minghelli et al. 2021; Zhu et al. 2020), and physical models must be verified by the actual water depth (Brando et al. 2009; Jay and Guillaume 2014). Due to data limitations, many studies have been limited to islands, lakes, and rivers. Therefore, passive water depth inversion was still limited by the lack of labeled data before the launch of the Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) in 2018, making it difficult to verify the experimental results (Neumann et al. 2019). The performance of ICESat-2 in shallow water sounding has opened up new possibilities for obtaining large-scale shallow water depth. Active and passive fusion sounding has become a new area of active research (Albright and Glennie 2020; Hsu et al. 2021; Liu et al. 2021b; Ma et al. 2020), which considers the water depth extracted by a single point as the true value, uses the empirical model to fit the model parameters, and then extends it to the relevant sea areas.

The premise of active and passive fusion sounding is to effectively filter ICESat-2 ATL03 data. The emitted energy of the photon-counting lidar system is much lower than that of the linear detection lidar, and the energy of the signal photons emitted and received by the system is relatively small and easily affected by various noises (Neumann et al. 2019). Several signal detection methods have been proposed for photon-counting LIDAR data (Chen et al. 2021a), including area-based methods (Awadallah, Abbott, and Ghannam 2014; Chen and Pang 2015), local statistics (Herzfeld et al. 2013; Magruder et al. 2012) and density clustering (Wang, Glennie, and Pan 2018; Zhang and Kerekes 2014). Since its public application, ICESat-2 has been widely used in various geographical areas including rivers and lakes (Armon et al. 2020; Brêda et al. 2019), ice sheets (Brunt et al. 2021; Brunt, Neumann, and Smith 2019), sea level (Buzzanga et al. 2021), shallow water depth (Hsu et al. 2021; Liu et al. 2021a; Ma et al. 2020; Parrish et al. 2019; Xu et al. 2021), land altimetry (Ye et al. 2021), etc. DBSCAN-based and OPTICS-based method is widely used for photons denoising in land surface. Adaptive density-based models can detect the ground and vegetation canopy in photon-counting laser altimeter data more accurately than DBSCAN in low-noise environments (Zhang and Kerekes 2014). There are also related studies on ATL03 denoising using this method (Gleason et al. 2021; Le Quilleuc et al. 2022; Wang et al. 2022), which relies on predefined parameters because of the high underwater noise and different underwater photon distributions that cannot be widely applied to different regions. Many studies have achieved promising results by improving DBSCAN (Liu et al. 2021; Ma et al. 2020; Xie et al. 2021; Xu et al. 2020). The main approach is to improve the empirical parameters, such as MinPts or Eps. However, these approaches are not universally applicable and may not produce new empirical parameters that are superior to the existing ones. Unlike DBSCAN, OPTICS-based does not explicitly create clusters but still relies on the input parameters (Tang et al. 2021; Zhu et al. 2020). Different methods have variable effectiveness for heterogeneous underlying surfaces and environments (Kui et al. 2023), such as the local statistics denoising method for reservoir surface(Xie et al. 2022), and spatial clustering and bimodal reconstruction method for forests (Huang et al. 2023), adaptive clustering and kernel density estimation for sea ice surface (Liu et al. 2023), self-adaptive denoising based on a genetic algorithm for land surface (Zhang et al. 2022), multilevel filtering algorithm for data with high background noise (You et al. 2023). ATL03 currently adopts a histogram-based filtering algorithm for products officially released by ATLAS. However, its performance is poor in areas with complex terrain areas (Neumann et al. 2021). Additionally, denoising methods exist for different surface types and observation conditions (Lao et al. 2023). However, the seabed signal weakens with increasing water depth and is affected by seabed reflectivity (Giribabu et al. 2023; Parrish et al. 2019). Therefore, denoising methods developed for land surfaces are not suitable for seabed photons extraction.

For the underwater topographic photon (UTP) extraction of ATL03 in shallow water, the photon distribution density on the seabed changes with changes in the water quality, seabed sediment, and seawater depth. The photon density varies significantly under different conditions, and the degree of dispersion varies with the change in seabed topography, which is unlike single-photon point cloud extraction on land. UTP extraction can be manually removed using an inefficient software platform (Babbel, Parrish, and Magruder 2021; Hedley, Velázquez-Ochoa, and Enríquez 2021). At present, the ATL03 UTP extraction method is primarily based on local statistics and density clustering, including AVEBM (Chen et al. 2021a), AEDTA (Wang et al. 2023), Density-based (Cao et al. 2021; Cao et al. 2023; Xie et al. 2023; Zheng et al. 2023), Neural Network-based (Meng et al. 2022) and DBSCAN-based method (Zhong et al. 2023). The AVEBM determines the filtering parameters based on the photon density distribution in different water environments and depths. It detects the signal photons on the water surface and underwater target point (UTP) (Chen et al. 2021a). The AVEBM has been applied in many studies (Chen et al. 2021b; Xie et al. 2021; Zhang et al. 2022a; Zhang et al. 2022b). However, AVEBM is suitable for the condition that the density of photons decreases linearly with water depth, which is not applicable to different seabed environmental conditions. AEDTA repeatedly calculates the height difference histogram of underwater photons and adaptively determines the threshold. However, it easily produces excessive denoising, particularly for weak beam data. These methods introduce empirical parameters or produce poor results, which are often determined manually and are difficult to apply to different environmental data (Wang et al. 2023). The density-based model is the most widely used algorithm for UTP extraction, such as density and distance-based method (Zheng et al. 2023), adaptive underwater point denoising algorithm considering the search direction and search size (Cao et al. 2023), improved local distance statistics method (Xie et al. 2023), and adaptive Gauss filtering technology based on the density of point cloud (Cao et al. 2021), these algorithms ignored characteristics of continuous underwater terrain. In addition, UTP extraction algorithms also include a method combining the DBSCAN algorithm and a two-dimensional window filter (Zhong et al. 2023), and a denoising algorithm based on BP neural network (Meng et al. 2022), which are not applicable in places with low signal-to-noise ratios (SNR). The differences in photon density in different regions must be considered to propose a generally applicable point-cloud denoising method.

The main difficulties in noise filtering underwater photon-counting data are as follows: First, it is difficult to obtain accurate results in areas with large slopes and inconsistent noise levels. Traditional filtering methods are challenging to apply to a wide range of data, and only achieve good results through manual adjustment of empirical parameters at the local level. This is due to several factors, including sunlight noise, variations in beam energy, and differences in water and sediment quality. Additionally, point-based filtering can easily filter sparse signals. Second, although ATL03 is widely used for UTP, there is still a lack of an evaluation standard reference for the denoising data of photon-counting technology sounding, and related research has only a confidence reference in the denoising process without a unified standard and framework (Ranndal et al. 2021).

This study proposes a general method for UTP extraction, which is an active contour method combined with a variable convolution kernel (ACVCK). This is the first underwater photon-filtering method that considers the adjacent photon density. To denoise seabed photons under different environmental conditions and photon densities, we propose an active contour with adaptive parameters that rasterizes subsurface photons. The optimal rasterization parameters were obtained using a kernel ridge regression (KRR) model and the photons were fitted to a curve as a whole. The corresponding parameter range was formulated according to the fitting results, and the confidence index of the ATL03 shallow-water sounding was marked. The innovation of the method includes a general method for underwater topographic photons (UTP) extraction, a new framework of photons denoising for underwater topographic profile, the profile signals extraction and evaluation method based on fitting idea, and a developing open-source UTP data set. This study aimed to propose a fully automatic production mode for ATL03 shallow-water data products that does not require manual adjustment of empirical parameters and provides a standard for the future application of photon-counting LiDAR data.

2. Research data and areas

2.1. Research areas

The UTP detected by the ATL03 data is mainly in a first-class water body, and the area within 40 m and the reflectivity of the seabed sediment are still acceptable. The data displayed in this manuscript are mainly from Robert Island (Figure 1c), Seven Connected Islets (Figure 1f), Oahu Island (Figure 1a), Discovery Reef (Figure 1e), Vuladdore Reef (Figure 1d), Yongxing Island (Figure 1b) in the South China Sea, and coastal areas (Figure 1g, h), with seabed signals in these areas, as shown in Figure 1, the source of the image is Sentinel 2A. Figure 1c shows the west of Yongle Atoll, which is located west of the Xisha Islands. Oahu is located northwest of the Hawaii Islands. Oahu is the third largest island among the Hawaii Islands of the United States and is located in the middle of the Pacific Ocean. Figure 1b, located in the Xisha District, Sansha City, Hainan Province, is the largest and most densely populated island in the South region. Figure 1f is located in the northeast of the Xuande and Xisha Islands and is the overall name of the large reef plate where the North Island and other islands are located. Figure 1d and e show the Yongle and Xisha islands, respectively. When the reef in Figure 1e ebbs, the entire reef can be exposed to the sea, and the reef in Figure 1d is an atoll.

Figure 1. The research areas and ATL03 data distribution involved in the paper.

2.2. Research data

(1) ICESat-2 ATL03

ICESat-2, launched in September 2018 and equipped with the Advanced Topographic Laser Altimeter System (ATLAS), is designed to measure the elevation changes of ice sheets and glaciers, the freeboard of sea ice, and the height of forests. It has also shown applicability for shallow water depth (Parrish et al. 2019). ATLAS uses a green (532 nm) laser and single-photon-sensitive detection (Neumann et al. 2019). The laser can operate at a high repetition rate of 10 kHz, which causes the laser emission interval along the track direction to be 70 cm (McGill et al. 2013). A diffractive optical element divides a laser beam into six beams, comprising three pairs of strong and weak beams, and one beam with a stronger pulse energy, approximately four times that of the other weak beams. The distance between each pair of beams is 90 m, and the distance between each pair of pairs of beams is 3 km. This arrangement provides greater spatial coverage (Markus et al. 2017). ATL03 contains the common noise of a single-photon-counting lidar system, mainly solar noise and atmospheric scattering (Herzfeld et al. 2013). In the ATL03 dataset (https://search.earthdata.nasa.gov/), a ‘confidence’ parameter (from 0 to 4) is used to classify each photon as a possible signal or noise. Errors such as atmospheric delay, earth tides, and system pointing deviation were corrected using the ATL03 dataset (Neumann et al. 2019). The time of ATL03 data used in this research is from 2018 to 2022, shown in Table 1 and Table 7, which include specific file names.

(2)	In-situ data of water depth

Table 1. Data of ATL03.

No.	Filename	D/N	Beam	Start	End	SNR	Kd(490) /m^−1
1	ATL03_20200419053723_03620701_005_01.h5	D	1l	16.169	16.182	0.44	0.094
2	ATL03_20210214031837_07961007_005_01.h5	D	2l	16.931	16.940	0.37	0.040
3	ATL03_20200419053723_03620701_005_01.h5	D	2r	16.162	16.177	1.30	0.010
4	ATL03_20210214031837_07961007_005_01.h5	D	1r	16.810	16.838	0.85	0.041
5	ATL03_20210617212611_12991107_005_01.h5	N	1l	16.321	16.347	5.74	0.036
6	ATL03_20210818183003_08571207_005_01.h5	N	2l	16.165	16.207	8.30	0.079
7	ATL03_20210716200216_03541207_005_01.h5	N	3r	16.329	16.358	15.62	0.043
8	ATL03_20210818183003_08571207_005_01.h5	N	2r	16.165	16.207	8.39	0.079

The experimental data included two types of water-depth data. The land and underwater topographic map of North Island in Figure 1f, 1:2000 scale measured in April 2013 (Ai et al. 2020). Measured data from Robert Island were collected in 2016. In this survey, we adopted an advanced full-waveform airborne radar system, using a 70 kHz pulsed Nd: YAG laser head to generate a green (532 nm) beam with twice the frequency, with an error of less than 0.2 meters (Su et al. 2018). The verification data for Oahu Island were collected using an airborne LiDAR SHOALS 3000 (Liu et al. 2021b). Figure 2 shows the data, the source of the image is Sentinel 2A. Because the positions of the in situ points do not completely correspond to the position of the ICESat-2 photons, the grid image was generated by kriging interpolation of the in situ point, and the pixel values of the grid image were extracted according to the latitude and longitude of the ICESat-2 ATL03 photons for comparison. In fact, it is challenging to perfectly match the measured data with ICESat-2 data. Several studies have used measured data interpolation methods to confirm this. Research has shown that the theoretical sounding accuracy of ICESat-2 is between 0.5 and 1 m. The experimental results of Section 4.4 show that the data in the three areas were all between 0.5 and 1 m. The interpolation results for these three areas are relatively reliable. Only North Island had a larger error, which was due to the coefficients of the sounding points in this area. Of course, the accuracy of the actual interpolation is not only related to the density of the measured data points but also to the degree of actual terrain relief. The interpolation error in this area was estimated to be within 1 m.

Figure 2. Example of the spatial distribution of in situ data.

3. Method

This paper proposes a general method combining the energy functional and nonlinear fitting methods for the UTP extraction of photon-counting LiDAR data and summarizes a new framework for such tasks aimed at producing a global UTP dataset. The technical framework of the study is shown in Figure 3 and mainly includes classification, signal extraction, confidence classification, evaluation, and application. First, a Gaussian statistical model was used to separate the photons from the sea surface and subsurface signals. The photons of the sea surface and above were removed, and all the photons of the ocean subsurface were obtained. Second, the parameter-adaptive active contour combined with a variable convolution kernel (ACVCK) method extracts the bottom signal contour and further extracts the underwater photons by rasterizing the subsurface photons and considering the pixels of adjacent photons. Third, the KRR was applied to fit the photon curve proposed by ACVCK and to calculate the optimal convolution kernel parameters by calculating the underwater signal density per unit area. Fourth, different buffer ranges were determined based on the fitting curve, and the underwater photon confidence index was marked in combination with the ACVCK extraction photons. Finally, we corrected the refraction and tidal errors of the filtered underwater photons to accurately solve the coordinates and compared the results with bathymetry data to verify the extraction accuracy and confidence index. We developed a website to display additional underwater photon data that can be downloaded (http://www.oceanread.com:5600/ExPress/ExDataPresentation).

Figure 3. Technical framework.

3.1. ACVCK for photon extraction

Currently, there is much research on UTP extraction methods. Point-based methods, such as elliptic and spherical filtering, are widely used. These methods are based on constructing a basic range to determine a photon label (Zhang and Kerekes 2014; Zhu et al. 2020). However, the water quality and bottom reflectance are quite different, and these parameters cannot adapt to different SNR, limiting their practical applications. The single-photon density is complicated for a detectable sea bottom, and there are difficulties in extracting sea bottom single photons, such as the density of signals being small, signals being limited by the detection depth, the number of reflected photons being different, and underwater after-pulse signals being strong. Most of these methods are either limited in their applicability or require extensive manual parameter tuning. Therefore, we propose the development of an automatic extraction method that is generally applicable to different water qualities and bottom reflectance. If this method can be applied and promoted, the continuous characteristics of underwater photons and the differences in the density of adjacent photons must be considered to mark the photon labels.

There are no manual empirical parameters in the ACVCK; therefore, the method was applied to different datasets. The photons were converted into raster by dividing the data horizontally and vertically. The signal photons were identified using the extracted range based on the energy difference. This method adaptively calculates the size parameters of the convolution kernel in ACVCK based on the KRR fitting model to obtain the maximum number of photons per unit length in the dead-time buffer of the fitting curve. Finally, we quantitatively estimate the signal extraction accuracy using data under different environmental conditions and propose an online open-source dataset. Figure 4 shows the entire process using data no.4 as an example. (1) $F_{\epsilon }(I_x,I_y,K_c)\stackrel{p_b,l_f}{⟶}max({\displaystyle \frac{n_{pb}}{le{n}_{lf}})\to I_x^{\mathrm{\prime }},I_y^{\mathrm{\prime }}\to F_{\epsilon }}$(1) (2) $n_{pb}=A_{pb}\cap B_{pb},A_{pb}=(x,y)\in in(C),B_{pb}=(x,y)\in b_f=l_f\pm di{s}_b$(2) (3) $le{n}_{lf}={\sum }_{i=0}^{n-1}\sqrt{{(x_{i+1}-x_i)}^2+{(y_{i+1}-y_i)}^2},n\in l_f$(3) (4) $di{s}_b=t\times C_w$(4)

Figure 4. Method process.

Formula (1) shows the core concept of ACVCK, $F_{\epsilon }$ is a function that combines the convolution kernel and energy functional parameters. The image is determined according to the horizontal axis distance ($I_x$) and vertical axis distance ($I_y$), where $K_c$ is the convolution kernel. Then, the seabed photons ($p_b$) are extracted using certain parameters, and the continuous profile $l_f$ is fit according to $p_b$. The length of $l_f$ is calculated as follows: $le{n}_{lf}$ shown in Formula (3), n represents all the points from $l_f$, and $x_i,y_i$ represents the horizontal and vertical coordinates of the photons. We calculate the intersection ($n_{pb}$) of the generated photon set ($A_{pb}$) using $F_{\epsilon }$ and the photons set in the buffer ($B_{pb}$) from $l_f$, as shown in Formula (2), where C is the fitting curve profile and $b_f$ is the buffer zone. Equation (4(4) $di{s}_b=t\times C_w$(4) ) calculates the size of the buffer ($di{s}_b$) based on the dead time of ATLAS (t) and transmission speed ($C_w$) of light in water. We determine the optimal horizontal and vertical axis distances when obtaining the maximum value of the ratio of $n_{pb}$ and $le{n}_{lf}$, that is, the number of extracted photons is the maximum for a unit length in the buffer zone when the optimal parameters are $I_x^{\mathrm{\prime }},I_y^{\mathrm{\prime }}$. Finally, the final seabed photon set is obtained according to the optimal parameters.

1. Separation for sea surface and sub-surface

The range of underwater photons must be determined before denoising the underwater photons, that is, separating the sea-surface photons. The sea surface height is calculated by the maximum value of photon density aggregation when the data is divided by 0.1 m in the vertical direction. A Gaussian model was constructed assuming that the photons conform to a Gaussian distribution within a range of 5 m above and below the sea surface, as shown in Formula (5), where h is the photon elevation within this range (Wang et al. 2023). To reduce the influence of photons on the sea surface, the photons on the sea surface and subsurface were separated using four times standard deviation ($\sigma$) to reduce the influence of high-density photons close to the sea surface, as shown in Formula (6), where $h_{bs}$ is the separation elevation between the sea surface and subsurface, and photons whose elevation is lower than $h_{bs}$ participate in denoising. Figure 4a shows the range of underwater photons involved in denoising, the results are marked in red. (5) $h=a{e}^{-{\displaystyle \frac{{(x-\mu )}^2}{2\times {\sigma }^2}}}$(5) (6) $h_{bs}=\mu -4\times \sigma$(6)

(2) Photon rasterization processing

Ocean subsurface photons are rasterized by calculating the number of photons within a certain range, including two parameters: horizontal and vertical intervals. Revising the distribution characteristics of photons in sparse areas is challenging when the statistical range is either too small or too large. A small range introduces noise into a unit, while a large range makes it difficult to determine whether the photons are seabed signals. The relevant parameters are rasterized according to a basic unit that determines the range of the convolution kernel. The number is mapped to the grey value range (0–255) after it is converted into a grid, as shown in Formula (7). $ras$ is the grey value of the grid; $N_{max},N_{min}$ are the maximum and minimum of the photon number within the grid range, respectively; $G{r}_{max},G{r}_{min}$ are the maximum and minimum of the grey range, which are 255 and 0, respectively. N is the number of photons per pixel. Figure 4b shows the rasterization results. (7) $\mathrm{ras}(\mathrm{x},\mathrm{y})={\displaystyle \frac{N_{max}-N_{min}}{G{r}_{max}-G{r}_{min}}\times (N-N_{min})}$(7)

(3)	Active Contour Generation

The Chan-Vese (CV) model is an active contour model (ACM), that was originally used for image edge detection and detects objects that define boundaries without gradients (Chan and Vese 2001). The ACVCK improves the CV model by providing new applications. This method minimizes energy by constructing an energy function, $F_{\epsilon }$, as shown in (8). Photon extraction is considered a special case of the minimum-partition problem. The application becomes a similar mean curvature flow in the level set, and the active contour containing the target photons stops at the desired boundary. However, it is difficult to achieve good results in seabed photon extraction due to the high noise and varying SNR. Therefore, a convolution kernel was introduced to improve the denoising accuracy. Since everything is interconnected, but nearby things are more closely related (Siewert et al. 2021), it is necessary to consider the influence of adjacent photons and a $3\times 3$ convolution kernel to represent the indirect performance of the adjacent threshold photons on the seabed photons, which enhances the signal and reduces the noise contribution. (8) $\begin{array}{rl}& F_{\epsilon }(I_x,I_y,K_c)=\mu {\int }_{\mathrm{\Omega }}\delta (\varphi (I_x,I_y)|\mathrm{\nabla }\varphi (I_x,I_y)|dxdy+\\ & {\lambda }_1{\int }_{\mathrm{\Omega }}{|R_{x,y}-c_{in}|}^{2}H(\varphi (I_x,I_y)dxdy)+\\ & {\lambda }_2{\int }_{\mathrm{\Omega }}{|R_{x,y}-c_{out}|}^2(1-H(\varphi (I_x,I_y)))dxdy\end{array}$(8) (9) $H={\displaystyle \frac{1}{2}(1+{\displaystyle \frac{2}{\pi }\arctan ({\displaystyle \frac{z}{\epsilon }))}}}$(9) (10) $R(x,y)={\displaystyle \frac{ras(x,y)\times K_c}{{\sum }_{i,j=0,0}^{i,j=d,d}{K}_{c_{ij}}}}$(10) (11) $K_c=\left|\begin{array}{ccc}a& & a\\ & a& \\ a& & a\end{array}\right|+\left|\begin{array}{ccc}& a& \\ a& a& a\\ & a& \end{array}\right|=\left|\begin{array}{ccc}a& a& a\\ a& 2a& a\\ a& a& a\end{array}\right|$(11)

Fixed parameter ${\lambda }_1={\lambda }_2=\epsilon =1$ in (8). $R_{x,y}$ was calculated using the convolution kernel $K_c$ (Equations 10(10) $R(x,y)={\displaystyle \frac{ras(x,y)\times K_c}{{\sum }_{i,j=0,0}^{i,j=d,d}{K}_{c_{ij}}}}$(10) and 1(1) $F_{\epsilon }(I_x,I_y,K_c)\stackrel{p_b,l_f}{⟶}max({\displaystyle \frac{n_{pb}}{le{n}_{lf}})\to I_x^{\mathrm{\prime }},I_y^{\mathrm{\prime }}\to F_{\epsilon }}$(1) 1). which is the pixel value after the convolution calculation. $c_{in},c_{out}$ are the average pixel values inside and outside the contour, respectively, using the level set function $\varphi$ to represent the curve. The length parameter μ for scaling is the number of columns along the track direction. H is the regularization function in Equation (9(9) $H={\displaystyle \frac{1}{2}(1+{\displaystyle \frac{2}{\pi }\arctan ({\displaystyle \frac{z}{\epsilon }))}}}$(9) ), where z is the energy value before the next iteration, which is used to calculate the Euler-Lagrange equation of $\varphi$. The method applies the ACVCK model to rasterized photon-counting LiDAR data, and the contour stops at the boundary between the signal and noise without considering the SNR, thereby obtaining the seabed signal. Figure 4c and d show the range of the seabed signal before and after rasterization.

3.2. KRR for parameter optimization

A seabed topographic profile exhibits a continuous feature that must be represented by a nonlinear function. KRR is a machine learning method with excellent nonlinear fitting capabilities. It combines ridge regression with kernel methods (Murphy 2012) to learn a linear function in the kernel-induced feature space of the data. This corresponds to a nonlinear function in the original input space for nonlinear kernels. The KRR computes the loss function using the squared error loss combined with l2 regularization. KRR is advantageous in situations where a nonlinear fit is required, or where the number of attributes exceeds the number of training examples.

All subsurface photon data cases were replaced with eigenvectors, that is $x_i\to v_i=v(x_i)$, where the number of dimensions may be much higher than the number of data cases. It can be interpreted that the algorithm is linear in the feature space; therefore, even if the dimensions of the feature space are much larger than the number of data cases, we do not need to access potentially infinitely long eigenvectors. In this application, the fitting value is the predicted value y of the test point x; $p_{bx},p_{by}$ are the coordinates of the seabed photons along the track and perpendicular to the track, respectively, and are calculated by w (Equations 12(12) $\mathrm{w}=(\lambda {\mathrm{I}}_d+v{v}^T{)}^{-1}vy=v(v^{T}v+\lambda {\mathrm{I}}_n{)}^{-1}y$(12) and 1(1) $F_{\epsilon }(I_x,I_y,K_c)\stackrel{p_b,l_f}{⟶}max({\displaystyle \frac{n_{pb}}{le{n}_{lf}})\to I_x^{\mathrm{\prime }},I_y^{\mathrm{\prime }}\to F_{\epsilon }}$(1) 3), where $\lambda$ controls the amount of shrinkage, and the larger the value, the larger the shrinkage, and the more robust the coefficient of collinearity. where K is the kernel function used in the fitting. The method uses the Laplacian kernel function, and $\gamma$ is the coefficient of the kernel function, as shown in Eq. (14(14) $\kappa (x)=K(x_i,x)=e^{-\gamma \Vert x_i-x\Vert },\gamma >0$(14) ). (12) $\mathrm{w}=(\lambda {\mathrm{I}}_d+v{v}^T{)}^{-1}vy=v(v^{T}v+\lambda {\mathrm{I}}_n{)}^{-1}y$(12) (13) $\begin{array}{rl}& y=krr(K,p_b,\lambda ,\gamma )\\ & ={\mathrm{w}}^{T}v(p_{bx})\\ & =p_{by}(v^{T}v+\lambda {\mathrm{I}}_n{)}^{-1}{v}^{T}v(p_{bx})\\ & =p_{by}(K+\lambda {\mathrm{I}}_n{)}^{-1}\kappa (p_{bx])}\end{array}$(13) (14) $\kappa (x)=K(x_i,x)=e^{-\gamma \Vert x_i-x\Vert },\gamma >0$(14)

3.3. Grading of UTP

UTP grading was performed using the results extracted by ACWCK and the buffer from the curve fitted by KRR. All subsurface photons were recorded as S, the result extracted by ACVCK was A, and KRR was used to construct a buffer related to the dead time, where the records in the double buffer were set as B₁, and the records in the double buffer were recorded as B₂, as shown in Formula (15). We divided the seabed photons into three levels: high, medium, and low. The photons are marked as high in both A and B₁. The photons marked as medians include two parts: one in B₁ not in A, and the second in A and between B₁ and B₂. The photons marked as low included two parts: one was between B₁ and B₂, but not in A, and the second was the photons in A but not in B₂. In summary, photons that are not in A but in the buffer drop the confidence to mark. Figure 4e shows the photons grading results, and the photons marked (a–e) correspond to the five situations of Formula (15), a-e represent the photon from $A\cap B_1$, $A\cap {\complement }_{B_2}{B}_1$, $B_1\cap {\complement }_{S}A$, $A\cap {\complement }_S{B}_2$, and ${\complement }_{S}A\cap {\complement }_{B_2}{B}_1$. (15) $\{\begin{array}{c}high,pts\in A\cap B_1\\ median,pts\in A\cap {\complement }_{B_2}{B}_1\cup B_1\cap {\complement }_{S}A\\ low,pts\in A\cap {\complement }_S{B}_2\cup {\complement }_{S}A\cap {\complement }_{B_2}{B}_1\end{array}$(15)

3.4. Photons signal post-processing

The extraction results from ACVCK were verified with in situ data, and post-processing was required to reduce errors. The ICESat-2 data photons were corrected for inherent errors using refraction correction and data acquisition time errors using tidal correction. The planar component of the refraction correction can be ignored when a first-order approximation is sufficient. Therefore, the corrected elevation for seafloor photon returns can be approximated using Equation (16(16) $z^{\mathrm{\prime }}\approx z+0.25416D$(16) ) (Parrish et al. 2019). In addition, tidal corrections were made to the in situ bathymetry data and ATL03 photons using the OTPS2 tidal model (Ma et al. 2020). After obtaining the corrected seabed photons, the KRR algorithm was used to fit the obtained photons, thereby obtaining a continuous topographic curve of the photons, and the difference between the sea surface height (see 3.1(1)) and the elevation value of the continuous topographic curve was used as the water depth. Figure 4f shows the corrected results. (16) $z^{\mathrm{\prime }}\approx z+0.25416D$(16)

4. Results

4.1. Qualitative validation

Typical data were selected to verify the accuracy and generalization of the ACVCK compared with existing methods, which should cover as many different SNR and environmental conditions as possible. Thus, we selected four types of data: D-L, D-R, N-L, and N-R. Two examples of each type of data are shown in this study, whose locations are displayed in Figure 1 and specific information in Table 1. More results are available in the open system we developed (the website is shown in Section 5.3). D/N indicates daytime and nighttime and Start/End refers to the latitude range. The manually labeled data for the presence of sea subsurface signals are shown in Figure 5 The subsurface photons were labeled as signals or noise with an SNR ranging from 0.3–16. In summary, the selected areas were representative. The results of the manual annotation are shown in Figure 5, where the parts are magnified to better visualize the seafloor topography through photon distribution (shown in the third column in Figure 5).

Figure 5. Label of UTP.

Classical signal photon extraction methods, such as the DBSCAN and OPTICS denoising methods (Zhang and Kerekes 2014; Zhu et al. 2020), are performed by constructing regions and counting the number of photons in the region. However, these methods are not suitable for underwater photons (Figure 8) and require uncontrolled empirical values during denoising, including the size of the range and number of photons in the region. Similar UTP extraction methods for sea subsurface signals, such as AVEBM (Chen et al. 2021a), construct elliptical filters that vary with depth, resulting in distance threshold parameters and parameter decay. However, the number of underwater photon signals does not strictly follow the distribution, which is delayed with water depth, limiting the performance of these methods. The relatively new method, AEDTA (Wang et al. 2023), performs UTP extraction based on the distance threshold between photons, which leads to excessive denoising, as shown in Figure 7, where the signals are incorrectly identified as noise, resulting in misclassification. Table 2 lists the empirical parameters.

Table 2. Empirical parameters.

No.	AVEBM	DBSCAN	OPTICS
1	0.01	1.5, 3	20
2	0.005	2.5, 3	15
3	0.02	1, 3	25
4	0.1	2.5, 3	15
5	0.04	3, 3	7
6	0.05	3, 3	10
7	0.2∼10	1, 3	15
8	0.1	3, 3	25

These are point-based denoising methods and the results are shown in Figure 8. It is apparent that such methods no longer satisfy the application requirements for sea subsurface signal photon extraction. The ACM has been widely adopted in image processing studies, mainly in two categories: edge- and region-based (Xie 2009). Due to the distribution characteristics of photons, which tend to cluster within a specific range, we utilized region-based methods. The results are shown in Figure 6.

4.2. Quantitative validation

The comparison in Figures 6–8 shows that ACVCK performs better than the existing methods with strong robustness. The denoising results of several methods were evaluated quantitatively using two parameters, F1 and IOU, as shown in Equation (17(17) $\{\begin{array}{c}F1={\displaystyle \frac{2P\times R}{P+R}\left(P={\displaystyle \frac{TP}{T_{ex}}R={\displaystyle \frac{TP}{T_{tr}}}}\right)}\\ IOU={\displaystyle \frac{P\times R}{P+R-P\times R}}\\ RMSE=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^N{(f(x_i)-y_i)}^2}{N}}}\\ MAE={\displaystyle \frac{{\sum }_{i=1}^N|f(x_i)-y_i|}{N}}\end{array}$(17) ). As shown in Table 3, point-based denoising methods have an advantage in filtering out photons that were originally parts of the signals, and ACVCK surpassed the other methods in terms of the two parameters above. (17) $\{\begin{array}{c}F1={\displaystyle \frac{2P\times R}{P+R}\left(P={\displaystyle \frac{TP}{T_{ex}}R={\displaystyle \frac{TP}{T_{tr}}}}\right)}\\ IOU={\displaystyle \frac{P\times R}{P+R-P\times R}}\\ RMSE=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^N{(f(x_i)-y_i)}^2}{N}}}\\ MAE={\displaystyle \frac{{\sum }_{i=1}^N|f(x_i)-y_i|}{N}}\end{array}$(17) Table 4 shows the data processing time (the unit is second), all the time is about several seconds, there is no obvious difference in efficiency for application. In fact, efficiency is related to many aspects, hardware performance, code program logic implementation, multi-threading, and multi-process all affect the efficiency of the algorithm. The advantages of ACVCK are automation and robustness.

Figure 6. Results of ACVCK.

Figure 7. Results of AEDTA.

Figure 8. Results of other methods.

Table 3. Comparison of methods.

No.	ACVCK		AVEBM		AEDTA		DBSCAN		OPTICS
	F1	IOU	F1	IOU	F1	IOU	F1	IOU	F1	IOU
1	0.89	0.81	0.57	0.40	0.61	0.44	0.43	0.28	0.44	0.28
2	0.84	0.72	0.42	0.26	0.45	0.29	0.49	0.32	0.59	0.42
3	0.89	0.80	0.75	0.60	0.69	0.53	0.82	0.70	0.61	0.44
4	0.88	0.78	0.58	0.41	0.60	0.43	0.70	0.54	0.65	0.48
5	0.92	0.86	0.82	0.70	0.85	0.75	0.61	0.44	0.77	0.63
6	0.93	0.87	0.77	0.63	0.86	0.75	0.66	0.49	0.77	0.63
7	0.92	0.86	0.89	0.80	0.88	0.78	0.92	0.86	0.57	0.40
8	0.92	0.86	0.89	0.80	0.86	0.75	0.91	0.84	0.72	0.56

Table 4. Efficiency comparison.

No.	ACVCK	AVEBM	AEDTA	DBSCAN	OPTICS
1	2	2	3	2	2
2	2	2	2	1	1
3	3	3	3	2	3
4	3	2	4	3	3
5	2	2	3	2	2
6	3	3	3	3	2
7	5	7	6	4	4
8	4	4	4	3	3

The extraction results were determined based on the successful retention of photons. However, our objective was to obtain a continuous seabed profile (seabed topography) through single-photon results; therefore, the extraction results were adopted as a quantitative evaluation from another perspective by comparing the results after photon fitting. Because of the difficulties in obtaining a stable and continuous fit using DBSCAN, AVEBM, and OPTICS, we fitted ACVCK and AEDTA and compared the results with the manual extraction results to calculate the MAE, Max, Med, and RMSE, as shown in Table 5. The extraction results of ACVCK were superior to those of AEDTA. It should be noted that the results retained most of the seabed photons compared with manual extraction. Therefore, the errors were mitigated, and the data were presented with high accuracy.

Table 5. Comparison of the fitting method.

No.	ACVCK				AEDTA
	RMSE	Max	Med	MAE	RMSE	Max	Med	MAE
1	0.25	1.37	0.07	0.14	0.48	2.99	0.10	0.26
2	0.66	2.65	0.09	0.31	0.43	0.91	0.33	0.37
3	0.14	0.87	0.03	0.07	0.52	3.44	0.07	0.21
4	0.34	1.73	0.13	0.22	1.46	5.14	0.53	0.95
5	0.39	2.73	<0.01	0.09	0.70	3.72	0.15	0.38
6	0.24	2.30	<0.01	0.06	0.54	4.04	0.12	0.27
7	0.28	3.06	<0.01	0.06	0.62	4.77	0.05	0.21
8	0.23	3.07	<0.01	0.05	0.30	3.45	0.04	0.12

The proposed method in this study relies on the significant grayscale difference between the signal and noise grids. However, we also increased the grayscale difference through the convolution kernel, followed by the fitting of the underwater terrain signal into an irregular curve. Therefore, the results are under the combined effect of these three conditions. The extraction accuracy of data with different SNR is generally higher for higher SNR than for lower SNR, such as day and night data. However, this relationship is not absolute. Additionally, the accuracy evaluation exhibited some fluctuations, which is also normal and to some extent expected, given that only eight data points were displayed. When the boundary between the signal and noise grids is blurred, if there are signals on the left and right sides, continuous signals can be extracted through fitting, which is also an advantage of this algorithm. However, if there are signals on the left and right, and the signal boundary is still blurry after convolution, it proves that there is indeed no signal in this area or that it is difficult to judge manually. Generally, as long as there is a signal, there can be differences through convolution amplification, and then, combined with fitting, the signal area can be extracted, and the signal photons can be determined. Generally, if humans can make clear judgments, they can obtain better results. In contrast, if humans cannot clearly discern whether a signal exists, it is possible that the photons have not penetrated the medium, which would make it difficult for the signal algorithm to make a judgment.

4.3. Elevation for grading

This study proposes a method for extracting UTP evaluations. To verify its reliability, we compared the accuracy (ratio of the number of correct photons at a confidence level to the total number of photons at this confidence level) of the three confidence levels with the manual extraction results, which can truly reflect the photon reliability at different confidence levels. The table indicates that the accuracy rate increases with increasing confidence level with a positive correlation, which demonstrates that the evaluation method is reliable. In Figure 6 and Figure 4e of the examples of the underwater noise classification results, we can deduce that the higher the degree of photon aggregation, the higher the degree of photon confidence, as shown in Table 6, where Hr, Mr, and Nr are the right photons in the label, and Ha, Ma, and La are the photons extracted by ACVCK.

Table 6. Reliability verification of evaluation methods.

No.	High		Median		Low
	Ratio	Hr/Ha	Ratio	Mr /Ma	Ratio	Nr /La
1	0.94	509/542	0.83	95/114	0.33	16/48
2	0.90	270/301	0.63	22/35	0.67	10/15
3	0.92	2060/2241	0.51	123/242	0.34	35/104
4	0.91	1190/1313	0.92	292/319	0.83	137/165
5	0.92	959/1039	0.67	24/36	0.65	17/26
6	0.93	2016/2159	0.79	27/34	0.69	9/13
7	0.93	6808/7299	0.87	317/364	0.68	19/28
8	0.93	6210/6647	0.82	318/389	0.73	96/132
Ave	0.92	–	0.76		0.62	-

4.4. Validation of in-situ data

ACVCK performed better than some existing methods in terms of denoising, compared with the results shown in section 3 .1. In this section, we post-processed the photons, including refraction and tidal corrections, and performed a correlation analysis using onboard high-precision underwater signals. All the data used are presented in Table 7 (No. 13 and No14 are for Section 5.3). Meanwhile, we achieved data interpolation and compared them with the ATL03 seabed results, owing to the limited resolution of the in situ data. Figure 9 shows the comparison between the ICESat-2 ATL03 seabed signal and the in situ bathymetry results. The first column shows the signal extraction results, the second is the comparison between corrected true values with ICESat-2 results, and the third is the correlation between the two types of data. ACVCK and AEDTA were comparable because they adopted the same validation data, whereas N, R2, and RMSE outperformed those of AEDTA (Wang et al. 2023).

Figure 9. Comparison with in-situ data.

Table 7. Comparison ATL03 for measured data.

No.	Filename	D/N	R/L	Start	End	Kd(490)/m^−1
9	ATL03_20181209231549_11050101_005_01.h5	D	2r	21.274	21.295	0.032
10	ATL03_20220304145339_11051401_005_01.h5	N	1l	21.28	21.30	0.032
11	ATL03_20220902061330_11051601_005_01.h5	N	3r	21.2851	21.312	0.072
12	ATL03_20221008163833_02741707_005_01.h5	D	3r	21.286	21.294	0.068
13	ATL03_20220317202009_13071401_005_01.h5	N	2l	16.51	16.5139	0.040
14	ATL03_20210214031837_07961007_005_01.h5	D	3r	16.9595	16.9625	0.041
15	ATL03_20181014192647_02470106_005_01.h5	N	3r	36.752	36.766	0.350
16	ATL03_20191028140111_04840501_005_01.h5	N	3r	21.435	21.445	0.276

The data time affects the accuracy of the results in some river estuaries or areas with greater human influence. However, our study area was located in the ocean, was less affected by human influence, and had no river influence. In addition, it may be affected by global sea levels; however, this effect is very small. The data in the three areas were all between 0.5 and 1, which is close to the theoretical accuracy.

5. Discussion

5.1. A new framework for UTP extraction

This study proposes a new pattern for photon extraction by fitting a plus-energy model. It transforms point cloud data into a raster map, converts the signal extraction problem into an image processing problem by extracting the signal range by combining the ACM with a variable convolution kernel, and obtains fitting results using the KRR. ACVCK can be generalized into a new generic framework for signal extraction, as shown in Figure 10, which demonstrates the process framework and code logic. The KRR algorithm plays a significant role in this framework by assisting in (1) determining the optimal range parameters, (2) determining the confidence level, and (3) fitting the continuous terrain results. For a continuous section of seafloor topography, the signals from one part of the seabed may be unclear, but those from adjacent areas may be opposite. Therefore, ACVCK has a clear advantage. It has no changing empirical parameters compared with traditional methods (DBSCAN and OPTICS) and is more stable with no excessive denoising compared with AEDTA.

Figure 10. Framework of photons for UTP. (a) Framework. (b) Logic of the methods.

5.2. Parameters calculation

1. Parameters

Other parameters are also included in this method. It is worth emphasizing that such parameters are fixed such that it is not necessary to change them as the data changes, which means that they are available for all data. During the ACVCK iterations, the step size was set to 0.2. The convolution kernel in ACVCK is used to fully account for the effects of domain photons, where a is 1 in Eq. (11(11) $K_c=\left|\begin{array}{ccc}a& & a\\ & a& \\ a& & a\end{array}\right|+\left|\begin{array}{ccc}& a& \\ a& a& a\\ & a& \end{array}\right|=\left|\begin{array}{ccc}a& a& a\\ a& 2a& a\\ a& a& a\end{array}\right|$(11) ).

The principle of ACM must be solved iteratively. During the solution process, regions with large differences in point cloud density can be difficult to handle consistently. Otherwise, low-density signal photons may be filtered out by high-density signal photons during the fitting process. Therefore, we divided the photon regions along the horizontal axis and calculated the number of photons within the sea surface range. To divide the raster, the minimum division unit along the track was set to three times the footprint size. We roughly estimated the density of single photons, where the highest photon density to the lowest photon density was less than twice the uniform processing range, and we performed iterative operations for each part separately to avoid excessive denoising caused by signal density differences. The division results of Data No.4 are listed in Figure 11, and it can be seen that statistical histograms with similar densities are divided into contiguous regions.

Figure 11. Area division.

We traversed the value of the statistical histogram (from left to right), the histogram value was added in an array one by one, the ratio of the maximum value to the minimum value in the array is calculated when a value is added to the array, until the ratio of the maximum value to the minimum value is greater than 2. Then the values in the array except the latest added value were one group, that is, the part between the two blue lines in Figure 11, which reduces excessive denoising of photon signals with different densities in the iterative process.

(2)

Iteration

Optimization and control of the model were achieved through iterative simulations of the energy model. To bring the function to convergence, we used a specific number of iterations for the energy-fitting model. However, too many iterations can waste computational resources, while too few iterations may not converge to the optimal solution. The convergence is judged by a parameter related to the length of the variable line ($\mu {\int }_{\mathrm{\Omega }}\delta (\varphi (I_x,I_y)|\mathrm{\nabla }\varphi (I_x,I_y)|dxdy$). For example, in Data No.4, divided into five regions (a-e), the decay in each region is shown. The change in the length parameter was obtained as the number of iterations increased. The number of iterations was determined when the derivative of the curve approached zero and was terminated. As shown in Figure 12, the degree of data aggregation and convergence number differed because of the different ranges of the five regions. Therefore, there are two ways to terminate the iterations: using a larger number of iterations, such as 500, but this would waste computational resources; therefore, we terminate the iterations when the fitted curve derivative tends to zero.

(3)	Effects of convolution range

Figure 12. Influence of iteration times.

We made a qualitative comparison by not adding kernels (that is 1 × 1 kernel) (Figure 13a), adding 3 × 3 kernel (Figure 13b), and adding 5 × 5 kernel (Figure 13c). Figure 13 shows that it is easy to miss signals without adding a kernel, as shown in a₁ and a₂; when the kernel is 5 × 5, the noise is easily introduced in areas with large terrain undulations, as shown in c₁; it will take the neighborhood photons into consideration when adding a kernel, the 3 × 3 convolutional kernel will also introduce a small amount of noise, as shown in b₁ and b₂, but it can be removed by optimizing the range parameter. The range parameters are 1/3 to control variables.

Figure 13. Influence of convolution kernel size.

The core of this method is the convolution range, which is determined by the number of photons per unit length of the fitted terrain obtained from KRR. The convolution range parameters must be obtained in two dimensions: the distance along the track and in the vertical direction. First, we consider the parameters associated with ICESat-2 in both directions, as shown in Eq. (18(18) $\{\begin{array}{c}{d}_z=c\times d_t\\ d_x=n_d\times d_{fp}\end{array}\Rightarrow \{\begin{array}{c}bi{n}_z=d_z\times C_{\mathrm{z}}\\ bi{n}_x=d_x\times C_x\end{array},C=0.1,0.2,\cdots ,1$(18) ). The dead time d_t for the vertical track of the Atlas data is 0.32 ns (Neumann et al. 2019) and the dead time distance is calculated, with c still taken as 3*10⁸ m/s and $d_z$ as 0.96 m as the ATL03 native data is not corrected for c in water. In the along-track direction, the influence of adjacent terrain is worth considering as the terrain is continuous; therefore, the convolution kernel dimension $n_d$ (set as 3 in this study) is taken into account, and $d_{fp}$ is 17 m. For data with different characteristics, specific values need to be determined for the different data, which are calculated by the maximum ratio of $n_{pb}$ and $le{n}_{lf}$, denoted as n/l. To control for variables, the two parameters of the KRR fit were set to 0.01.

The effects of different $bi{n}_z$ and $bi{n}_x$ on n/l are shown in Figure 14, with the horizontal and vertical axes, C, ranging from 0.1–1. The boxed grid has a better value (PV) of n/l. It can be observed that the better values of $C_{\mathrm{z}}$ are mainly distributed between 0.3 and 0.9, and the better values of $C_x$ are mainly distributed between 0.1 and 0.5. As in Equation (19(19) ${\displaystyle \frac{d_{fp}}{n_d}<bi{n}_x<d_{fp}\times n_d\times 0.5,{\displaystyle \frac{d_z}{n_d}<bi{n}_x<d_z}}$(19) ), the ranges of $bi{n}_z$ and $bi{n}_x$ should contain at least 1 d after convolution and should not exceed 1.5 d. Finally, the optimal C is determined according to the parameter range and the maximum ratio of $n_{pb}$ and $le{n}_{lf}$, so that the final result is generated. (18) $\{\begin{array}{c}{d}_z=c\times d_t\\ d_x=n_d\times d_{fp}\end{array}\Rightarrow \{\begin{array}{c}bi{n}_z=d_z\times C_{\mathrm{z}}\\ bi{n}_x=d_x\times C_x\end{array},C=0.1,0.2,\cdots ,1$(18) (19) ${\displaystyle \frac{d_{fp}}{n_d}<bi{n}_x<d_{fp}\times n_d\times 0.5,{\displaystyle \frac{d_z}{n_d}<bi{n}_x<d_z}}$(19)

(4)	Parameter fitting

Figure 14. Influence of convolution kernel range.

A nonlinear fitting algorithm was selected to fit the seabed single-signal photons due to the complex underwater topography and irregular undulating conditions. Figure 15 illustrates the effects of the different methods on the fitting results. The chi2 kernel function fits lines that are too smooth to represent continuous profiles in great detail by observing the fitting results for data Nos. 1, 4, 5, and 8 under different environmental conditions. The rbf kernel function is prone to large errors in small areas of missing data; therefore, the Laplace kernel was selected for KRR fitting. The Laplace kernel performs well for this task, as shown in the figure. Figure 16 shows that the two variables have a significant effect on the fitting results. Therefore, it is essential to determine the best-fitting parameters. We obtained the optimal parameters using the GridSearchCV auto-tuning method (Pedregosa et al. 2011), where $\lambda \in ({10}^{-4},{10}^0),\gamma \in ({10}^{-4},{10}^4)$, and at 10-fold intervals.

Figure 15. Fitting method.

Figure 16. Fitting parameters.

5.3. Method robustness

We added the UTP extraction results under Class-II water body data (No 15 and 16 in Table 7 and Figure 1), shown in Figure 17. The parameters used to evaluate water quality are Kd(490) data from NASA Ocean Color (https://oceancolor.gsfc.nasa.gov/l3/). The Kd(490) range for all experimental data is approximately between 0.032∼0.35, which is shown in Table 1 and Table 7. In fact, water quality will affect the ultimate detection depth of photons. As long as it is a detectable signal, it will leave obvious signal photons, the algorithm will achieve good results.

Figure 17. Results of data with different K_d.

ACVCK achieves extraction of seabed single-signal photons under different environmental conditions (day and night, strong and weak beams) and SNR. Numerous experiments were conducted to verify the robustness of this method. Besides the results shown in this study, we have developed a website (http://www.oceanread.com:5600/ExPress/ExDataPresentation/) to show more data (the data results are being expanded). The aim was to form an open dataset for seabed single-signal photon data and make it available for download by relevant research scholars. The website will also be gradually improved and more features will be developed in the future, including autonomous data uploading and data visualization. In addition, we shared the results and data of this study for review via GitHub (https://github.com/iwzhcode/ACVCK), and the source code will be made publicly available at this address.

5.4. Improvements and directions

In this study, we propose a new framework for single-signal photon extraction from a seabed that has room for improvement. (1) We iterated the energy model separately by zoning, which was subject to chance errors and may still have intrazone energy differences, leading to excessive denoising in some areas. A superior ACM model should be established, where the continuous terrain is regarded as a whole and local small-scale energy differences are considered to achieve better results. (2) When some regions are too noisy, the results can be optimized by improving the low-confidence range, such as data Nos. 9 and 10 in Table 7, by adjusting the low criterion $A\cap {\complement }_S{B}_2$ to $A\cap {\complement }_{B_3}{B}_2$, thereby removing regions farther away from the buffer. Currently, a subsea single-signal photon extraction method is under development.

Engineering applications of single-photon bathymetry remain challenging. Single-photon bathymetry is used in engineering applications and includes the following steps: (1) Obtaining the surface height during surface signal extraction; (2) obtaining the surface height below the surface; (3) acquiring the area where the seabed signal exists; (4) removing noise to extract the signal; and (5) correcting the seabed signal. Currently, there is limited research on step (3), and studies have shown that carrying other payloads can solve this issue. However, different payloads have varying detection depth limits, and therefore, further investigation is required. Currently, ICESAT-2 is a 6-beam system that can only detect depths of up to 40 m, and the development of more beams and higher-energy single-photon lidars onboard is the next step.

6. Conclusion

This study proposes a new extraction method for UTP images that combines active contours with a variable convolution kernel. The proposed method overcomes the limitations of traditional methods, such as indeterminate empirical parameters, poor adaptability, and susceptibility to excessive denoising. It is a general mathematical model framework that can be applied to a variety of denoising problems. The initially extracted underwater photons were fitted using the KRR, and the optimal parameters were determined using the maximum number of photons per unit length to ensure the universality of the method. An innovative method for evaluating UTP extraction was proposed using the fitted curve buffer and regional photon results to classify the signal as high, medium, or low. This method achieved excellent performance for data in different environmental situations, including day/night and strong/weak data. The F1 of the different data was above 0.84 compared with the manually labeled data, the IOU was above 0.72, and the RMSE of the fitted data was basically within 0.5. The correlation coefficients (R2) were all above 0.94 compared to the in situ data. We eventually made the data and method open-source, generated a small batch of underwater signal data, and made it publicly available to illustrate the stability of the method. This dataset is still under expansion, and our goal is to produce a complete UTP dataset to facilitate the development of shallow-water bathymetry.

Acknowledgements

We are grateful to the NSIDC for providing the ATL03 data, and thanks to ATL03 data for their contribution to the related study area. We also would like to thank the editors and reviewers for their valuable comments on the paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available in an online website, at http://www.oceanread.com:5600/ExPress/ExDataPresentation. The results supporting the paper are available at https://github.com/iwzhcode/ACVCK, and the code will be available at this site.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China (grant numbers, 41930535, 62071279, 41871382, 42301501 and 42001416), and the Independent Research Program of Key Laboratory of Land Satellite Remote Sensing Application, MNR (grant number BN2302-6).