ICESAT-2 Shallow Bathymetric Mapping Based on a Size and Direction Adaptive Filtering Algorithm

The US Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) satellite adopts 532 nm single-photon lidar with shallow sea bathymetry capability. In order to realize high-precision and automated shallow sea bathymetric mapping based on ICESat-2 photon data, an adaptive underwater point denoising algorithm that considers the search direction and search size is proposed in this article, and a detailed data processing process is discussed to further verify the technical feasibility. The accuracy of direct bathymetry and active–passive fusion bathymetry from ICESat-2 is systematically analyzed using airborne in situ data. First, the underwater photon points are separated by surface position identification; then, the signal point cloud search strategy is improved, the search size increases with water depth, the search angle is rotated and the direction of the maximum number of point clouds is taken as the main direction, and the threshold is automatically determined by histogram Gaussian fitting of the point cloud density to achieve automatic underwater signal extraction; then, the refraction correction is carried out based on the light geometry and the water depth is obtained; finally, the active–passive fusion bathymetry is performed by combining the optical remote sensing images WorldView-2 and Sentinel-2, and the accuracy is verified by using the airborne lidar bathymetric data provided by NOAA. The experimental results show that the proposed denoising algorithm can accurately discriminate the underwater signal/noise, and the overall accuracy is better than 86%; the root-mean-square error (RMSE) of ICESat-2 direct bathymetry is between 0.42 and 0.98 m; the RMSE of active–passive fusion bathymetry is between 0.84 and 1.88 m. Our workflow and experimental results demonstrate a means of using ICESat-2 to produce relatively accurate bathymetric maps in shallow, clear water environments.

marine environmental protection, marine military operation, and other occasions [2], [3], [4]. Remote sensing technology has a high acquisition efficiency and good data timeliness when compared to traditional sonar bathymetry, making it an alternative means in shallow sea bathymetry.
Remote sensing bathymetry technology can be active lidar or passive optical imagery. Airborne-based and satellite-based are the most common acquisition platforms. The most common ones are airborne LIDAR bathymetry (ALB) and satellite optical image bathymetry. ALB has developed rapidly in recent years, and the most advanced CZMIL system can obtain integrated point clouds of water and land with the accuracy of (3.5 + 0.05d) m (2σ, σ is the standard deviation, and d is water depth) in plane and (0.302 + (0.013d) 2 ) 1/2 m (2σ) in elevation [5], meeting the demand for high-precision coastal zone mapping. However, the ALB system is expensive, and only a small number of users can afford it, limiting the application scope. The satellite-derived bathymetry includes empirical models [6], [7], [8], physical models [9], [10], [11], and two-media photogrammetry [12], [13], [14]. The empirical model relies on in situ depth data to solve model parameters, and it cannot be used in areas where depth control information is unavailable; the physical model requires accurate water body parameters and bathymetric accuracy is closely related to water body parameter accuracy; and the two-media photogrammetry requires a calm water surface, or it will destroy the three-dimensional ray geometry structure. In conclusion, while remote sensing images are large in scope and easy to obtain, bathymetric accuracy is affected by a variety of environmental factors, and the relative error of bathymetry is frequently greater than 10%, which can only meet general water depth investigations.
Single-photon lidar is a new type of lidar developed in recent years. It primarily uses single-photon detectors as receiving devices, such as avalanche photodiodes operating in Geiger mode and photomultiplier tubes. These detectors are two to three orders of magnitude more sensitive than traditional linear detection lidar, and they make it easier to achieve direct 3-D imaging with micropulses, high refrequency, and multiple beams [15], [16]. The Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) launched in 2018 uses single-photon detection mode [17], [18], and the still-planned Lidar Surface Topography satellite also intends to use single-photon lidar [19].
Shallow water bathymetry is not the main engineering goal of ICESat-2, but it has bathymetric capabilities due to the use of a green laser in the 532 nm band [20]. On the one hand, ICESat-2 can be used as a direct bathymetric measurement of shallow waters with a maximum penetration depth of about 1 Secchi (about 40 m in clear water) and a root-mean-square error (RMSE) of better than 0.6 m in bathymetry [21]. On the other hand, laser points can be used as control points to perform fusion bathymetry with passive optical images and thus obtain regional bathymetric maps. The bathymetric accuracy is related to the fusion model and the RMSE can reach about 1 m [22], [23], [24], [25], [26]. ICESat-2 bathymetry mainly includes steps of underwater signal point identification, refraction correction, depth calculation, reference correction, and active-passive fusion. Refraction correction [23], [27] and active-passive fusion model receive more attention, and automatic identification of underwater signal point with high accuracy is much more challenging. ICESat-2 Algorithm Theoretical Basis Document [28], [29] provides Poisson denoising algorithm and differential, regressive, and Gaussian adaptive nearest neighbor filtering (DRAGANN) algorithm, which do not consider the photonic properties of the water body and cannot be directly applied to bathymetric processing. Ma et al. [30] used DBSCAN algorithm to extract underwater signals, Chen et al. [22] investigated variable-length search kernels, and Hsu et al. [31] used segmented median filtering for denoising, all of which are more effective in processing specific data but still rely on empirical threshold setting and require manual intervention to achieve satisfactory results [25].
In this article, we use ICESat-2 data in shallow sea area to carry out point cloud denoising, and use WorldView-2 (WV2) and Sentinel-2 (SA2) optical images for active-passive fusion bathymetry, and carry out accuracy verification based on the in situ data provided by NOAA. The main objectives include the following.
1) To provide a fully automatic high-precision underwater point denoising algorithm to get rid of human intervention. 2) To provide a more detailed photonic bathymetry data processing process to further verify the technical feasibility. 3) To carry out systematic accuracy evaluation of ICESat-2 direct bathymetry and active-passive fusion bathymetry, providing reference for applications. The rest of the article is organized as follows. Section II describes the theoretical approach, including water surface identification, size and directional adaptive filtering (SDAF), refraction correction, active-passive fusion model, test area, and data. Section III describes the description of the experimental results. Section IV discusses some specific issues. Finally, Section V concludes the article. Fig. 1 depicts the data processing flow of this article. First, the elevation histogram is counted in the vertical direction, and the water surface elevation is calculated by using Gaussian fitting method to divide the point cloud into above-surface points, surface points, and underwater points; subsequently, a noise-removal algorithm with adaptive search length and search direction is designed for underwater points to achieve accurate identification of signal points; then refraction correction is carried out based on the light geometry, and the mapping depth is obtained through datum correction; finally, the accuracy of ICESat-2 direct bathymetry and active-passive fusion bathymetry is evaluated in comparison to an airborne lidar data set.

A. Calculation of Water Surface Elevation
As shown in Fig. 2(a), the ICESat-2 data are profiled point clouds distributed along the track. The laser pulse passes through the air-water-seabed and is affected by environmental noise. Return signals may be recorded throughout the transmission process. For subsequent accurate identification of the underwater signal, as well as for calculating the water depth, the water surface elevation needs to be calculated first. Because the number of photons returned by the laser pulse at the water surface is much greater than the background noise and underwater signal, the frequency of photon occurrence in each elevation section on the point cloud profile can be counted to form an elevation histogram [as shown in the blue section of Fig. 2(b)]. The water surface location is determined by taking the elevation with the highest frequency of occurrence in the elevation histogram.
The real ocean environment has surges, tides, etc., resulting in the actual sea surface not being a perfect plane. The real ocean environment has surges, tides, etc., resulting in the actual sea surface not being a perfect plane. Due to the water penetration of 532 nm photons, many photon points are returned in a certain range above and below the water surface. These points are Gaussian distributed in the elevation direction. In order to identify the water surface points and water surface elevation more reliably, this article performs a Gaussian fit to the elevation histogram [as shown by the red dashed line in Fig. 2(b)]. The mean value of the Gaussian parameter is taken as the water surface height and  the points within ±3 times the variance marked as water surface points.

B. Size and Direction Adaptive Filtering
Due to the effect of background noise and dark currents, multiple echo signals may be generated from the same emitted signal, resulting in a low signal-to-noise ratio for single photon point clouds. In comparison to the noise, the localization of the ground echoes has a higher density. The vast majority of single-photon LIDAR filtering algorithms [32], [33], [34], [35] rely on local density as their primary filtering criterion. To do this, they search within the neighborhood of the target point p i and count the density value D i (expressed as the number of points), if D i > T (T is the threshold), it is a signal, otherwise it is noise. The main challenge is to set the proper neighborhood range for the distribution characteristics of underwater photons and to calculate the threshold automatically to ensure that the underwater signal can be extracted completely.
The main characteristics of underwater photons are shown in Fig. 3: 1) the number of photons decreases significantly with increasing water depth, with higher density at very shallow water (kernel 1) and only a few points at the deepest depth (kernel 3); 2) the underwater topography changes, and the main direction of photon distribution is not always horizontal. Inspired by the papers [36] and [37] and DRAGANN algorithm [29], this article proposes an SDAF algorithm, in which the search range increases with water depth, the rotating search direction takes the largest number of points as the main direction, the denoising threshold is automatically determined by Gaussian fitting, and the residual isolated points are removed by quadratic segment fitting. The specific steps of the algorithm are as follows.
1) Determine the search size. Since photon points have aggregation characteristics in a certain direction, define an elliptical search kernel centered at point p. The distance between any point q and point p is denoted as where x is the distance along the flight direction, corresponding to the horizontal axis of Fig. 3 (also see Figs. 7 and 8); h denotes the elevation, corresponding to the vertical axis of Fig. 3; a and b denote the long and short axes of the ellipse, respectively; x, h, and dist(p, q) are in meters; when dist(p, q) < 1, it means that point q is within the ellipse of point p.
The ellipse size are designed to vary with depth, a p = a 0 + 2.5dh, b p = 0.1a p , dh is the instantaneous water depth corresponding to point p, and a 0 = 10 m. The scaling parameters of the variable dimensions are set empirically, a p = 10 m, b p = 1 m at the shallowest point, a p = 110 m, b p = 11 m when dh = 40 m.
3) Determine the threshold by density histogram Gaussian fitting. Calculate the density of each point according to the size and direction adaptive ellipse, and then count the number of occurrences of each density to form a density histogram. The noise is located on the left side of the density histogram and the signal is located on the right side. Multiple Gaussian functions are used to fit the noise and signal waveforms to determine the signal-to-noise distinction threshold, and the relevant technical details can be found in papers [29], [38]. 4) Rejection of residual noise. After the aforementioned processing, the discovered signal points may still contain cluster outlier noise (see Fig. 4). We conduct segmented quadratic fitting in the along tract direction to further eliminate the noise. The quadratic function is where h is geodesic height; x represents along track distance; and c 0 , c 1 , and c 2 denote the quadratic function parameters. The specific steps are as follows. 1) Segment all the "signal points" to be processed in the direction of the track, the empirical value of the segment length in this article is 200 m. 2) Use all the "signal points" within the segment to fit the quadratic polynomial, and calculate the unknown parameters c 1 , c 2 , and c 3 . 3) Calculating the distance from the current point to the fitted curve point by point, if the distance is greater than the threshold T, the current point is marked as noise. 4) Looping steps 1)-3), stopping when no noise is identified or the number of loops is greater than 3 times, the empirical threshold T for each loop is 20, 10, and 5, respectively.

C. Refraction Correction
The ICESat-2 ATL03 level point cloud represents the location where the laser propagates in a straight line without considering underwater refraction. Parrish et al. [21] present a concise refraction correction method. Fig. 5 depicts the light underwater geometry, where the water surface is assumed to be an ideal plane. The horizontal deviation of the observation point A and the real point A due to refraction can be ignored for the satellite platform, and only the vertical deviation is corrected [12], [39].
In Fig. 5, n 1 and n 2 are the refractive indices of air and water, respectively, taking the values of 1.0 and 1.34. θ 1 and θ 2 are the angles of incidence and refraction, respectively, and the ground angle ref _elev is provided in the ATL03 file, θ 1 = π/2 − ref _elev. θ 2 is calculated based on the law of refraction n 2 /n 1 = sin θ 1 / sin θ 2 . The point A corresponds to the water depth z, which is obtained by making a difference between the elevation value of A and the water surface elevation (Section 1.1). In the triangle ST z, S = z cos θ 1 . The distance R is also obtained based on the law of refraction n 2 /n 1 = S/R, and the distance P and angle α are solved in triangle RPS using the cosine theorem and the sine theorem, respectively where φ = θ 1 − θ 2 . Since β = γ − α, the amount of refraction correction in the vertical direction can be calculated using Δz = P · sin β.

D. Active-Passive Fusion Bathymetry Model
Optical remote sensing bathymetry establishes the relationship between the electromagnetic radiation incident on an optical sensor for a discrete water sample unit and its corresponding water depth based on the radiative transfer equation or an empirically-derived equation, and uses the reflectance value for Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. each pixel of an optical image to calculate the water depth value at each corresponding water pixel. In this article, the Lyzenga multiwavelength inversion model [6] is used where z is the water depth; L(λ i ) is the off-water radiation corresponding to the band λ i ; L ∞ (λ i ) is the off-water radiation corresponding to the deep water area of the band λ i ; α i is the model parameter; and N is the number of bands used.
In order to carry out bathymetric inversion, absolute radiance conversion, atmospheric correction, sun-glint elimination [40], and water-land separation are also required for remote sensing images. Both ATLAS point clouds and optical remote sensing images are projected from the geodetic coordinates (longitude, latitude, and geodetic height) under the WGS84 ellipsoid to the UTM plane to realize the alignment of point clouds and image pixel. The water depth obtained from the ATLAS point cloud is taken as true value and substituted into (6) using the least squares method to calculate α i . Once α i have been determined, the offwater radiation of each pixel can be used to obtain pixel-by-pixel water depth for shallow areas.

E. Bathymetric Accuracy Evaluation
The bathymetry obtained by ICESat-2 point cloud is the instantaneous bathymetry relative to the sea level at the moment of imaging, and tidal correction is required for subsequent applications and accuracy assessment. If the acquisition time of ICESat-2 and satellite image are both known, in order to perform active-passive fusion bathymetry, the relative elevation difference between them can be calculated by querying the tide table or with the help of tidal model (such as NAO.99b model [41]). If the measured verification data have been converted to a chart bathymetric datum, the ICESat-2 bathymetry needs to be converted accordingly before comparing the accuracy.
For the accuracy assessment, a total of three metrics, RMSE, relative bathymetric error (RBE), and R 2 (R 2 = r 2 , where r is the correlation), are used in this article where h i represents the actual water depth at point i and h i represents the water depth measured by ICESat-2 at point i.
The smaller the RMSE, the greater the accuracy. The smaller the RBE, the greater the accuracy. However, when water depth becomes shallower, the smaller the error will cause a larger RBE, so we need to pay attention. R 2 indicates the proportion of the total sum of squares of deviation that can be explained by the regression square. The value is between 0 and 1, and the closer to 1 means the higher the accuracy.

F. Test Area and Data
As shown in Table I  ICESat-2 laser data and optical image specific information are listed in Table I. For example, St. Croix 20191015-2r indicates that the laser data was taken on October 15, 2019, and the ground track number of laser beam used is 2r. The latitude range was introduced to focus on the data processing near the coastal zone and reduce the amount of data. In terms of optical images, only Qilianyu and Yongle Atoll have WV2 images, whereas other regions have only SA2 images. The WV2 satellite is operated by the U.S. company MAXAR. WV2 image is a four-band multispectral data set with spectral bands of 450-510 nm (blue), 510-580 nm (green), 630-690 nm (red), 770-895 nm (near infrared), and a geometric resolution of 2 m. WV2 images have high geolocation accuracy with CE90 < 3.0 m (CE90 means circular error at 90 probability) without ground control points, and can be used directly without geometric correction. The SA2 satellite is developed and operated by the European Space Agency. SA2 has 13 spectral bands and different geometric resolutions; here only, B2 (blue, 458-523 nm), B3 (green, 543-578 nm), B4 (red, 650-680 nm), and B8 (NIR, 785-900) with 10 m resolution are used. SA2 geometric positioning accuracy can reach CE90 < 20 m, with a direct positioning error of no more than two pixels relative to its 10 m geometric resolution, and is also used directly without considering geometric correction.
For validation data, St. Croix, St. Thomas, Guam, Saipan, and Honolulu all use 1 m resolution digital elevation models (DEMs) provided by NOAA. The ALB methodology was applied in each of these areas to measure water depth. The bathymetric data in St. Croix and St. Thomas were acquired by RIGEL VQ88O-GII equipment in January 2013 with a vertical accuracy of ±18 cm [42]. In Guam and Saipan, the data were acquired by Hawkeye 4X equipment in February 2020 with a vertical accuracy of ±10 cm [43]. Honolulu data were gathered in February 2020 using CZMIL, with horizontal accuracy of (3.5 + 0.05d) m and vertical accuracy of (0.30 2 + (0.013d) 2 ) 1/2 m. SDE-28S+ single-beam sonar was used in May 2017 to measure depth of Qilianyu with a vertical accuracy of (1 ± 0.1%d) m. A total of 27 631 points were collected with a maximum water depth of -47 m. Yongle Atoll is very large (about 280 km 2 ), and only the Ganquan Island in its northwest corner has an actual water depth data measured by SHOALS-3000 in January 19, 2013. The horizontal and vertical accuracy are about ±2.5 m and ±15 cm, respectively. Due to the high point density, the 1m resolution DEM was also produced for later use.

A. Results of Adaptive Filtering
According to the methodology described in this study, surface elevation computation, underwater signal identification, and refraction correction were carried out. Fig. 7 displays some of the data processing results. In order to quantify the accuracy of the denoising algorithm, the signal points manually interpreted by three operators and cross-checked with each other were used as the real signals. The underwater signals, particularly the sparse signals in deep water, may be successfully extracted from the visual effect by SDAF. The automatically recognized signal locations are remarkably consistent with the manually interpreted signal places. For St. Thomas 20220818-3l, the water depth reaches about 45 m (before refraction correction) on the right side of 9.68 × 10 5 m in the horizontal axis, and the signal points are sparsely distributed. The SDAF's adaptive design allowed for the successful identification of this portion of the weak signal. Before refraction correction, the leftmost side (near 1.884 × 10 6 on the horizontal axis) of Qilianyu 20210413-3r data is about 46 m deep, and the SDAF successfully extracted the sparse weak signal.
The SDAF algorithm's extraction accuracy was quantitatively evaluated by comparing it to the DRAGANN [29] and adaptive variable ellipse filtering bathymetric method (AVEBM) [22] algorithms, and the results corresponding to seven data in Fig. 7 are shown in Table II. When the underwater photon aggregation is high, DRAGANN, a denoising algorithm developed for terrestrial scenes, can theoretically extract part of the underwater signal. The AVEBM algorithm also adjusts the search range based on the depth of the water, and dynamically sets the denoising threshold at different depths based on experience; whether the experience value is appropriate has a direct impact on denoising accuracy. Table II where A is the 3 × 3 matrix in Table II Croix 20191015-2r with 79%, 89%, and 91% overall accuracy, St. Thomas 20200818-3l with 77%, 85%, and 91%, and Saipan 20220323-3l with 62%, 76%, and 86%, respectively. The denoising accuracy of all three sets is DRAGANN<AVEBM<SDAF. The main reasons are the large water depths, the large differences in the distribution density of photons, and the more pronounced topographic undulations. These factors determine the denoising effect of "variable size + variable direction (SDAF)" > "variable size (AVERBM)" > "fixed direction+ fixed size (DRAGANN)." For Yongle Atoll 20190222-1l, all three algorithms agree in their accuracy of 97%. This is due to the shallow depth in this region and the high level of photon point aggregation. The two comparison algorithms DRAGANN and AVEBM agree with the results of [22] and [29]. Fig. 8 depicts a comparison of the three algorithms' results for Saipan 20220323-3l. DRAGANN can only extract dense signal points at shallow water depths. While AVEBM correctly identifies the sparse signal points to the right of 1.689 × 10 6 m, it fails to identify the inclined sparse points between 1.686-1.689 × 10 6 m.  For data Guam20210130-2r, AVEBM is the least effective, with an accuracy of only 67%. This is primarily due to the empirically determined dynamic threshold being insufficiently reasonable, which works well with other data but performs poorly with this data. DRAGANN performs better than SDAF, owing to the abundance of shallow water points in section 1.4700 − 1.4705 × 10 6 . SDAF misses some of the signal points in this segment due to the use of a quadratic fit. However, at the deepest position corresponding to 1.469 × 10 6 , DRAGANN is unable to extract the sparse signal points, whereas SDAF successfully identifies them.
The three algorithms show the same denoising ability in data Qilianyu 20210413-3r and Yongle Atoll 20190222-1l, with an accuracy of 92% and 97%. This is primarily due to the submerged substrate's high reflectivity in these two regions, as well as the high photon density and low topographic relief. It should be noted that although Fig. 7(f) and (g) appears to have more pronounced topographic undulations, this is actually owing to the visual effect caused by the range of the horizontal axis display. The underwater photon slope varies relatively little when shown 1:1 on the horizontal and vertical axes, allowing DRAGANN and AVEBM to produce successful extraction results.

B. Results of ICESat-2 Bathymetry
The direct bathymetric accuracy of the ICESat-2 photon point cloud was verified using ten sets of data from St. Croix, St. Thomas, Guam, Saipan, and Honolulu. The Qilianyu area had sparse data and the Yongle Atoll area did not have ICESat-2 groundlines passing directly through Ganquan Island, so the accuracy of these two areas could not be directly compared. Table III displays the direct bathymetric accuracy results, and Fig. 9 depicts the scatter plot of ICESat-2 measured bathymetry versus validated bathymetry. It should be noted that the signal points in Table III are the automatic extraction result of SDAF, and the number is not necessarily equal to Table II. For example, the number of St. Thomas20200818-3l is 301 in Table II and 138  in Table III. This is because the geographic range of the airborne verification data (see Table I) here is smaller than the extracted ICESat-2 point could, and the laser point outside the range needs to be removed. The green solid line in Fig. 9 indicates the 1:1 straight line, and the red solid line indicates the regression fitted straight line of ICESat-2 bathymetry and the validation bathymetry. The bathymetric RMSE in Table III ranges  The correlation coefficients R 2 of all ten sets of data exceed 0.962, which indicates the high accuracy of ICESat-2 bathymetry points, and it is generally considered that the two sets of data are extremely strongly correlated when 0.8 < r < 1.0. For RBE, the minimum is 4.7% (St. Croix20191015-2r) and the bigger reaches 23.2% (Guam20210228-1r) and 24.2 (Honolulu20190511-3l).
The different colors in the scatter plot of Fig. 9 respond to the different point densities. St. Croix20191015-2r has a total of 1210 points with a maximum water depth of -19.4 m. From Fig. 9(a), the distribution of point clouds in different water depths is more uniform, and the maximum point density occurs in -14   Guam20210228-1r has a total of 731 points with a maximum water depth of -13.1 m, which are basically distributed from 0 to -1 m [see Fig. 9(f)]. Honolulu20190511-3l has 365 points with a maximum water depth of -26.8 m, whereas most points are distributed from 0 to -3 m. For those two datasets, a smaller absolute error will produce a larger RBE, so the 23.2% and 24.2% in Table III are reasonable. In two sets of data, we successfully extracted signal points at -46 m (Saipan20220323-3l) and -45 m (Honolulu20190115-1l). Fig. 9(h) and (i) shows that, even though there are only a few points in the deepest water, they are all near the 1:1 regression line with no obvious offset. This demonstrates the effectiveness of the denoising algorithm in this article, which can extract weak signals in deep water.

C. Results of Active-Passive Fusion Bathymetry 1) Fusion Model Selection:
Active-passive fusion models have a large impact on bathymetric results, and although there are comparative analyses of fusion models in the literature, whether the relevant conclusions are applicable to linear control points like ICESat-2 needs further study. In this article, a comprehensive comparison of three common models is conducted under the same experimental conditions. Although the prediction ability of the model is related to the number of control points, it does not mean that the more control points, the higher the model accuracy. It is more important that the control point data are evenly distributed across depths to ensure the prediction accuracy of the model in different water depths. Usually, ICESAT-2 point clouds are denser in very shallow water (e.g., 0-3 m), and the point clouds become thinner when the depth becomes larger. In this article, before fusion bathymetry, the ICESAT-2 signal points are stratified by depth, and the dense point cloud in the very shallow water region is diluted, and all the signal points in the deeper region are retained; subsequently, one-half of the points are randomly selected as control, and the other half are used as check points, which can be seen that the number of check points in Table IV is often less than half of the signal points in the previous paper (see Table II).
The three tested models are the multiband model, the twoband ratio model, and the BP neural network (BPNN) model. Among them, the multiband model is shown in (6), and the two-band ratio model is referred to the literature [44]. In this article, the BPNN model uses a 3-layer network, and the number of nodes in the input layer is 3, which corresponds to the first three values of the six band ratios (b1/b2, b1/b3, b1/b4, b2/b3, b2/b4, b3/b4) of the image element after the principal component transformation, and the number of nodes in the output layer is 1, which corresponds to the water depth value of the image element. Table IV shows the internal compliance accuracies of the three models. The RMSE of the multiband model ranged from 0.634 to 1.904 m, R 2 ranged from 0.833 to 0.960, and the minimum RBE was 11.4%; the RMSE of the ratio model ranged from 0.605 to 2.137 m, R 2 ranged from 0.681 to 0.947, and the minimum RBE was 12.2%; the RMSE of the BPNN model ranged from 0.834 to 2.939 m, with R 2 ranging from 0.546 to 0.945 and a minimum RBE of 14.3%. Taking St. Croix 20191015-2r as an example, each accuracy metric showed that multiband model is better than the BP model, and the BP model is better than the two-band ratio model. Overall, the bathymetric accuracy of the multiband model was superior to that of the ratio model, with only two data sets (Honolulu 20190511-3l and Yongle Atoll 20190222-1l) having RM SE Ratio < RMSE M ultiband . Compared with the ratio model, BPNN model does not show obvious advantages. Among the 14 groups of data, seven groups of RM SE Ratio are smaller, and the other seven groups of RM SE BP N N are smaller. Fig. 10 depicts the bathymetric accuracy of each set of data as a bar graph, where the three widths, ranging from big to tiny, represent the multiband model, BP model, and dual-band ratio model, respectively, and the three colors, red, yellow, and green represent RMSE, R 2 , and RBE, respectively. From Fig. 10, the multiband model usually has smaller RMSE, smaller RBE, and  R 2 closer to 1. To investigate the reasons, only 2 bands of image data are used in the ratio model, the amount of input information is too small. The BPNN model has many influencing parameters, the number of nodes, the function type of nodes, etc., which will affect the bathymetric inversion performance of the model, and repeated experiments are needed to establish the best neural network. These are the main reasons for the poor accuracy of the dual-band model and the BP model compared with the multiband model.
The scatter plots of two datasets employing the three models are shown in Fig. 11 (St. Croix 20191015-2r and Saipan 20200924-1l). The check points for the multiband model [ Fig.  11(a) and (f)] are evenly distributed on both sides of the red fitted line, whereas those for the two-band model [ Fig. 11(b) and (e)] are obviously shifted to one side in the deep-water area, and those for the BP model [ Fig. 11(c) and (f)] are similar. This suggests that the multiband model performs better for the experimental data in this article than the two-band model.
Integrating the above experimental conclusions, the multiband model is selected below to produce the bathymetric map of the experimental area and carry out the subsequent analysis. The same control points are also used to create the bathymetric maps.
2) Comparison of Different Image Accuracy: The optical image band setting, resolution, atmospheric correction algorithm, and other factors will affect the active-passive fusion bathymetry results. To investigate this issue, we collected two types of images, SA2 and WV2, in two experimental areas of Qilianyu and Yongle Atoll (shown in Table I), using the same atmospheric correction algorithm, to observe the relationship between image resolution and bathymetry results. Table V shows the results of the comparative bathymetry experiments. The bathymetric accuracy does not improve with increasing spatial resolution, and the RMSE of SA2 and WV2 are nearly the same, though some SA2 data have smaller value (e.g., Qilianyu-1r, Yongle-3l). If we plot the scatter plots (not provide here), we can find that there is no significant difference between two type imagery for Qilianyu 20210413-3r. For Yongle Atoll 20190222-1l, WV2 image has more large error points deviating from the fitted straight line than SA2, but the overall accuracy is very close.
Considering data accessibility, free SA2 images should be preferred when carrying out multiregional large-scale bathymetry, and will not affect bathymetric accuracy due to resolution. For local areas, WV2 can provide bathymetric maps with higher resolution, which is more advantageous in engineering applications.
3) Absolute Accuracy of Fused Bathymetry: All ICESat-2 laser points in each test region were used as bathymetric control points. The absolute bathymetric accuracy of active-passive   Fig. 11, and the plots have been diluted due to the large number of points. Except for some obviously deviated points in the upper left side of Fig. 12(a), the points in the remaining images are evenly distributed on both sides of the fitted line, and the fitted line is relatively close to the 1:1 solid line in the three scatter plots. Fig. 13 depicts six active-passive fused bathymetry pseudocolor map created with SA2 images, demonstrating that the water depth is highly consistent with the subjective visual perception of multiband images.

A. SDAF Design Ideas and Usage Limitation
The SDAF method is designed to effectively collect deepwater weak signal points and establish the denoising threshold without the need for human input. The search size rises with the water depth and rotates the search direction in order to keep  the deep-water weak signal as much as feasible. The maximum number of points is then noted as the density corresponding to the target point.
As shown in Fig. 14(a), the SDAF algorithm judges the obvious noisy point p to be initially marked as a signal. This is due to the large water depth of point p, which corresponds to a large search range, and a total of 12 points located within the point p search ellipse after rotation, with density D p = 12. The density histogram generated by SDAF is shown in Fig. 14(d). Two Gaussian functions representing the signal (red line) and noise (green line) are used to optimally fit the histogram. The coordinate of horizontal axis corresponding to the intersection of two Gaussian functions is the denoising threshold T [T = 8.3 in Fig. 14(d)], and a detailed solution can be found in [29]. Since D p > T , the point p is judged as a signal. Such misjudgment leads to discrete noise in the initial SDAF processing results, which is mainly identified and removed by the quadratic topographic fitting in the algorithm.
By comparing Fig. 14(c) and (d), it can be found that there is a significant change in the SDAF density histogram compared with DRAGANN with fixed search length and direction. The deep-water sparse points after SDAF processing have larger density values and are more likely to be labeled as signals compared to DRAGANN. Fig. 14(c) also shows that DRAGANN may not be effective because the histogram does not necessarily show the "bimodal" characteristics like land data due to the high-density variation of underwater signal points. DRAGANN may be effective when the underwater points are densely distributed with little density variation (e.g., Qilianyu20210413-3r, Yongle Atoll 20190222-1l).

B. Active-Passive Fusion Bathymetry Limitation
Although ICESat-2 single-photon Lidar itself has high bathymetric accuracy, the fusion processing with optical images will inevitably introduce errors, and the success of fusion bathymetry is also closely related to the bathymetric distribution of ICESat-2 points, remote sensing waveband response, resolution, etc. [24], [45].
Using Saipan Island as an example (see Fig. 15), this article employs two sets of laser data: 20200924-1l is located on the right side of the image, primarily in the shallow water area above -10 m; 20220323-3l is located on the left side of the image, primarily in the deeper water area below -10 m; the airborne verification bathymetry data contains all the deep-water shallow water area. When the active-passive fusion model is built using the 2020924-1l shallow water point cloud, the model meets the requirement of internal conforming accuracy [see Fig. 15(b)]. Due to the high proportion of shallow water points and the lack of deep-water control points below -10 m, the model's generalization ability decreases sharply, and the predicted values at depths below -10 m deviate almost completely from the true values, as shown in Fig. 15(c).
The SA2 image in this region limits the use of the deepwater data 20220323-3l. The multiband bathymetric model (6) requires the specification of the deep water off-water radiance L ∞ (λ i ) and the satisfaction of L(λ i ) > L ∞ (λ i ). For NIR band, L ∞ (λ NIR ) is usually 0. This SA2 image has too many pixels with L(λ NIR ) = 0, leading to the failure of multiband bathymetric model construction. The above phenomenon may be related to the resolution of SA2 or the setting of NIR band range. This issue is not found on the WV2 images of Qilianyu and Yongle Atoll, both of which satisfy L(λ i ) > L ∞ (λ i ).

V. CONCLUSION
This article proposes an ICESat-2 single photon point cloud denoising algorithm SDAF based on adaptive search size and direction. It can automatically determine the search size with water depth and take the direction of maximum laser point density as the main direction. SDAF avoids the disadvantages of the traditional model parameter extraction method, which requires manual intervention, and can, to some extent, solve the sparse and weak signal point extraction problem in deep water. The test results on several sets of ICESat-2 data from six regions, such as St. Croix, St. Thomas and Guam, show that the overall denoising accuracy of the proposed SDAF algorithm is better than 86% compared to manually extracted bathymetry points. The two conventional algorithms used for comparison achieve only 62% and 76%. Our proposed method is better than the conventional algorithm in terms of denoising accuracy and deep-water point extraction effect.
In this article, ICESat-2 point clouds, optical remote sensing images, and bathymetric verification data from several survey areas are used to perform single-photon bathymetry accuracy verification and active-passive fusion shallow-sea mapping. After denoising, refraction correction, and tidal correction, the RMSE of ICESat-2 point cloud bathymetry is between 0.42 and 0.98 m. The fusion with optical images can obtain shallow sea depth pixel by pixel, with some loss in bathymetric accuracy, and the overall RMSE is between 0.84 and 1.88 m. The multiband model performs relatively well on the active-passive fusion model. In terms of fused bathymetric accuracy, the sentinel image with 10 m resolution is not worse than the WV2 with 2 m resolution. The next work needs to investigate the fusion bathymetric model with more adaptability and less error.