Retrieval of Nearshore Bathymetry around Ganquan Island from LiDAR Waveform and QuickBird Image

Featured Application: The work presented in this paper provides methods for the retrieval of nearshore bathymetry based on active and passive remote sensing, and will serve as a reference for shipping safety, coastal construction, island development and utilization, etc. Abstract: Optical remote sensing is an e ﬀ ective means of water depth measurement, but the current approach of mainstream bathymetric retrieval requires a large amount of onsite measurement data. Such data are hard to obtain from places where underwater terrains are complicated and unsteady, and from sea areas a ﬀ ected by issues with rights and conﬂicts of interest. In recent years, the emergence of airborne light detection and ranging (LiDAR) provided a new technical means for ﬁeld bathymetric survey. In this study, water depth inversion was carried out around an island far from the mainland by using remote sensing images and real LiDAR waveform data. Multi-Gaussian function ﬁtting was proposed to extract water depth data from waveform data, and bathymetric values were used as control and validation data of the active and passive combination of water depth inversion. Results show that the relative error was 5.6% for the bathymetric retrieval from LiDAR waveform data, and the accuracy meets the requirements of ocean bathymetry. The average relative error of water depth inversion based on active and passive remote sensing was less than 9%. The method used in this study can also reduce the use of LiDAR data and the cost, thus providing a new idea for future coastal engineering application and construction.


Introduction
Water depth is one of the basic elements of seabed topographic mapping and marine environment survey. The key and basic data provided by water depth monitoring are used in coastal engineering construction, ship navigation, and island development and utilization. The traditional bathymetric method of using ship-borne echo sounding achieved remarkable success. However, survey vessels and surveyors cannot easily carry out the traditional methods in shallow water (0-2 m) due to the influence of tide and submerged reef, indicating a restriction of these methods. Compared with traditional onsite acoustic measurements, the remote sensing inversion of water depth has the advantages of synchronous measurement in large coverage areas, high efficiency, and low cost. Remote sensing is an effective complement to traditional bathymetric methods in places where underwater terrains are complicated and unstable, and in sea areas affected by issues with rights and conflicts of interest.

Related Work
The bathymetric system of LiDAR matured as a technology in recent years. However, the data processing software of the bathymetric system is limited to some extent, and bathymetric accuracy is affected by the processing of waveform data, which are affected by bottom reflection, suspended matter, and sea waves [8]. The bathymetric algorithm of the LiDAR waveform accordingly became a research focus of many scholars. The bathymetric inversion methods mainly include the peak detection method [9,10], mathematical function fitting method [11][12][13], and convolutional filtering method [14][15][16]. Billard [17] proposed a method to detect seabed reflection signal by removing the low-frequency part of the echo signal through high-pass filtering. Wagner [9] employed the peak detection method to search the echo signal from the surface and the bottom of the sea, and the method was found to be feasible. Wong [18] decomposed two kinds of independent reflection signals, namely, sea surface echo signals and sea bottom echo signals, from which water depth was then estimated. Chen [19] improved the inversion accuracy of water depth by using the phase ambiguity resolution technique to correct wave and tide measurements. Yao [20] established a LiDAR echo waveform model based on the bi-Gaussian pulse fitting method and applied it to separate surface and bottom reflection signals. A minimum depth of 0.4 m was retrieved from their study. Zhang [21] utilized a Gaussian function entailing the double fitting of simulated waveform data to retrieve water depth. Their results showed that a mean absolute error (MAE) of 15.6 cm and a mean relative error (MRE) of 4.58% could be achieved at the water depth of 1-15 m using the LiDAR depth-detection model. In practice, a portion of the LiDAR effective signal leaches and the low-frequency signal is removed when using a filtering method, but the error is larger and the efficiency rate is slower when a peak detection method is used directly. The resultant fitting effect is unsatisfactory for deep waters due to the large energy gap between the echo signals from the sea surface and the sea bottom. Moreover, the bathymetric inversion error for deep waters is greater than that in shallow waters.
The accuracy of depth inversion by passive remote sensing is lower than that by LiDAR. Meanwhile, bathymetric retrieval by active optical remote sensing has high precision and maneuverability. Consequently, bathymetric retrieval methods based on active and passive remote sensing were studied in recent years. Pacheco [22] utilized the multi-band logarithmic linear model and Landsat 8 remote sensing image to invert shallow water depths (0-12 m) in areas around islands. The inversion algorithm of passive optical remote sensing was optimized using LiDAR bathymetric data. The average error of the inversion results was 0.2 m, and the median error was 0.1 m. Tian [23] carried out bathymetric retrieval by active and passive remote sensing with the use of a passive optical image from Landsat-8 and airborne LiDAR depth data, and subsequently analyzed the results of different density LiDAR data for depth inversion by passive optical sensing. Their results showed that the change in LiDAR density only minimally affected the water depth inversion results, and the average relative error changed by 0.3%. Many scholars also conducted substantial research and experiments on sediment classification and topographic mapping using active and passive remote sensing data [24][25][26]. Depth inversion by active and passive remote sensing can not only improve the accuracy of passive optical depth inversion, but also reduce the use of active optical data, as well as reduce costs.
On the basis of the above discussions, this study takes Ganquan Island as an example and uses the bathymetric values extracted from Aquarius LiDAR waveform data. Multi-Gaussian function fitting is performed to derive the control points and validation points, and the bathymetric method of combining active and passive remote sensing is carried out in the sea area around islands far from the mainland.

Research Area
Ganquan Island (111 • 35 10" east longitude and 16 • 30 28" north latitude) belongs to the Xisha Islands in China. The island has a tropical monsoon climate, and it is approximately 500 m from east to west (width) and 700 m from north to south (length), with an area of approximately 0.3 km 2 . The waters around Ganquan Island are clear and limpid and largely unaffected by human activities. The study area and multispectral image are shown in the Figure 1a,b respectively.
Appl. Sci. 2019, 9,4375 3 of 17 The accuracy of depth inversion by passive remote sensing is lower than that by LiDAR. Meanwhile, bathymetric retrieval by active optical remote sensing has high precision and maneuverability. Consequently, bathymetric retrieval methods based on active and passive remote sensing were studied in recent years. Pacheco [22] utilized the multi-band logarithmic linear model and Landsat 8 remote sensing image to invert shallow water depths (0-12 m) in areas around islands. The inversion algorithm of passive optical remote sensing was optimized using LiDAR bathymetric data. The average error of the inversion results was 0.2 m, and the median error was 0.1 m. Tian [23] carried out bathymetric retrieval by active and passive remote sensing with the use of a passive optical image from Landsat-8 and airborne LiDAR depth data, and subsequently analyzed the results of different density LiDAR data for depth inversion by passive optical sensing. Their results showed that the change in LiDAR density only minimally affected the water depth inversion results, and the average relative error changed by 0.3%. Many scholars also conducted substantial research and experiments on sediment classification and topographic mapping using active and passive remote sensing data [24][25][26]. Depth inversion by active and passive remote sensing can not only improve the accuracy of passive optical depth inversion, but also reduce the use of active optical data, as well as reduce costs.
On the basis of the above discussions, this study takes Ganquan Island as an example and uses the bathymetric values extracted from Aquarius LiDAR waveform data. Multi-Gaussian function fitting is performed to derive the control points and validation points, and the bathymetric method of combining active and passive remote sensing is carried out in the sea area around islands far from the mainland.

Research Area
Ganquan Island (111°35'10" east longitude and 16°30'28" north latitude) belongs to the Xisha Islands in China. The island has a tropical monsoon climate, and it is approximately 500 m from east to west (width) and 700 m from north to south (length), with an area of approximately 0.3 km 2 . The waters around Ganquan Island are clear and limpid and largely unaffected by human activities. The study area and multispectral image are shown in the Figure 1a and 1b respectively.

Aquarius LiDAR Waveform Data
The Aquarius LiDAR bathymetry system (Optech Company, Canada) is an amphibious measuring system for shallow sea areas which can measure both land information and shallow water information within 15 m. The bathymetry system uses a blue-green laser (532 nm). The specific parameters are shown in Table 1. The waveform data obtained by an Aquarius-airborne LiDAR were acquired from the findings of an onsite survey conducted on 19 January 2013. Flight time was approximately 3 h, and total data volume was approximately 55 GB. The information stored in each waveform data included reflected pulse signals, received pulse signals, and GPS time. The left side of the yellow vertical line in Figure 2 represents the transmitting pulse, while the right side represents the receiving pulse. Different targets or regions appear with varying return signal characteristics, mainly in terms of energy value, shape, and quantity of echo signals.

Aquarius LiDAR Waveform Data
The Aquarius LiDAR bathymetry system (Optech Company, Canada) is an amphibious measuring system for shallow sea areas which can measure both land information and shallow water information within 15 m. The bathymetry system uses a blue-green laser (532 nm). The specific parameters are shown in Table 1. The waveform data obtained by an Aquarius-airborne LiDAR were acquired from the findings of an onsite survey conducted on 19 January 2013. Flight time was approximately 3 h, and total data volume was approximately 55 GB. The information stored in each waveform data included reflected pulse signals, received pulse signals, and GPS time. The left side of the yellow vertical line in Figure  2 represents the transmitting pulse, while the right side represents the receiving pulse. Different targets or regions appear with varying return signal characteristics, mainly in terms of energy value, shape, and quantity of echo signals.

QuickBird Passive Optical Image
The QuickBird satellite provides a panchromatic band with 0.61-m spatial resolution and a multi-spectral image with 2.44-m spatial resolution, as shown in Table 2. The passive optical image adopted in this study was the multispectral image ( Figure 1b

QuickBird Passive Optical Image
The QuickBird satellite provides a panchromatic band with 0.61-m spatial resolution and a multi-spectral image with 2.44-m spatial resolution, as shown in Table 2. The passive optical image adopted in this study was the multispectral image ( Figure 1b  The auxiliary data used in this study mainly included tide tables for the period covering 2012 and 2013 published by National Marine Scientific Data Center. The time period of the LiDAR data measured onsite was approximately 3 h, and the tides changed considerably during those periods. Therefore, performing a tidal correction on each water depth value converted from the waveform data was necessary to ensure the accurate use of control points and validation points.

Data Processing and Method
For the LiDAR waveform data, we propose a multi-Gaussian function fitting method to extract water depth. The bathymetric values were converted into mean lower low water (MLLW) data to unify the standards of measurement. A series of preprocessing steps for the QuickBird image, such as radiation calibration and atmospheric correction, were needed. The improved logarithmic conversion ratio model (ILCRM) was used for the active and passive combination of water depth inversion in this study.

Radiation calibration
The ground object signals received by satellite sensors are recorded with a dimensionless Digital Number (DN) value, and the process of converting the DN value into a radiation brightness value with practical physical significance is called radiation calibration. The basic principle of radiometric calibration is to establish the quantitative relationship between the DN value and the radiation brightness value in the corresponding field of view.
The calibration formula for the remote sensing image of QuickBird is where L is the inputted pupil radiation brightness of the sensor, DN is the pixel gray value, absCalFactor is the absolute scaling factor, and effectiveBandwidth is the effective width of the spectrum. The information on absolute calibration factor and effective spectral width can be found in the calibration file (*.IMD). The details are listed in Table 3. Table 3. Radiation scaling transformation parameters of QuickBird.
Band absCalFactor effectiveBandwidth 1.542420 × 10 −2 1.14 × 10 −1 After radiometric calibration, the QuickBird remote sensing image is converted from one with the original dimensionless DN into one with a radiance value with µW/cm 2 × nm × sr units.

Atmospheric correction
Satellite sensors are affected by the scattering and absorption of atmospheric particles, aerosols, and molecules in the process of receiving solar photoelectric and electromagnetic radiation. The effects of scattering attenuate the signals received by the sensors. Meanwhile, some of the scattered light from the sky directly penetrates the satellite sensor without interacting with ground objects. This part of the signal does not contain any ground object information; we call this phenomenon atmospheric path radiation. Therefore, the atmospheric effects need to be eliminated to obtain highly realistic ground radiation or reflectance.
The atmospheric correction method used in this study was the 6S model method. The 6S model, which is the abbreviation of second simulation of the satellite signal in the solar spectrum, was developed by adding a new scattering algorithm based on the 5S model and, thus, can better reflect the real situation than the LOW resolution atmospheric TRANsmission (LOWTRAN) and MODerate resolution atmospheric TRANsmission (MODTRAN) models. Figure 3 shows the technical route of the 6S atmospheric correction that consists of three parts: Parametric setting of the atmospheric correction model, radiation calibration, and atmospheric correction [27].
Appl. Sci. 2019, 9,4375 6 of 17 The atmospheric correction method used in this study was the 6S model method. The 6S model, which is the abbreviation of second simulation of the satellite signal in the solar spectrum, was developed by adding a new scattering algorithm based on the 5S model and, thus, can better reflect the real situation than the LOW resolution atmospheric TRANsmission (LOWTRAN) and MODerate resolution atmospheric TRANsmission (MODTRAN) models. Figure 3 shows the technical route of the 6S atmospheric correction that consists of three parts: parametric setting of the atmospheric correction model, radiation calibration, and atmospheric correction [27]. The parameters of the 6S atmospheric correction model were set as follows: ground elevation, solar observation geometry, satellite observation geometry, and other parameters obtained directly from the *.IMD file corresponding to the metadata of the remote sensing image; the spectral response function obtained from the data disclosed by the sensor; the atmospheric model parameters with default values obtained according to the remote sensing image data, mainly in the aerosol mode, atmospheric mode, and solar spectrum function; and a customization of the actual situation. We selected marine aerosol and tropical atmospheric mode in accordance with the needs of the study area.
After radiation calibration, the DN value of the remote sensing image is converted into an apparent radiation brightness value, followed by a conversion into an apparent reflectance value by Equations (2) and (3). Apparent reflectance is the ratio of radiation brightness at the entrant pupil of the satellite sensor to the descending sun brightness on top of the atmosphere, which also includes the role of the atmosphere. The parameters of the 6S atmospheric correction model were set as follows: Ground elevation, solar observation geometry, satellite observation geometry, and other parameters obtained directly from the *.IMD file corresponding to the metadata of the remote sensing image; the spectral response function obtained from the data disclosed by the sensor; the atmospheric model parameters with default values obtained according to the remote sensing image data, mainly in the aerosol mode, atmospheric mode, and solar spectrum function; and a customization of the actual situation. We selected marine aerosol and tropical atmospheric mode in accordance with the needs of the study area.
After radiation calibration, the DN value of the remote sensing image is converted into an apparent radiation brightness value, followed by a conversion into an apparent reflectance value by Equations (2) and (3). Apparent reflectance is the ratio of radiation brightness at the entrant pupil of the satellite sensor to the descending sun brightness on top of the atmosphere, which also includes the role of the atmosphere.
where ρ TOA is the apparent emissivity, d is the relative distance between the sun and the earth, and ESUN λ is the average solar radiation value of the band with the unit W/m 2 . The value of ESUN λ can be calculated according to the spectral response function of the satellite sensor. The calculation formula is where λ 1 and λ 2 are the spectral wavelengths of the sensor with the unit µm, and E(λ) is the solar radiation value outside the atmosphere with the unit W/ m 2 µm . Equation (3) can be calculated from the solar spectrum response function provided by the World Radiation Center. After the model parametric setting is completed, the 6S atmospheric correction model can be run to obtain the atmospheric correction parameters of X a , X b , and X c . Subsequently, the apparent reflectance data of the image are converted into surface reflectance data. The conversion formulas are as follows: where ρ s is the surface reflectance, and X a , X b , and X c are the conversion parameters calculated by the 6S atmospheric correction model. The desired surface reflectance data were, thus, obtained.

Principle of bathymetry based on passive optical remote sensing
Sunlight has a perspective in a water column. The signal received by the optical remote sensor contains reflection information of the sea bottom to the sunlight, and this information is the physical basis of the passive optical remote sensing inversion of shallow water depth. In other words, the information reflected from the seabed to the optical remote sensor is the direct reflection of underwater terrain, that is, the information source of optical remote sensing bathymetry. The attenuation coefficient of sunlight in water determines the penetrable depth of the light in water. The smaller the attenuation coefficient of light in water is, the better its transparency in water is. This phenomenon is the basic principle of water depth inversion.
In the process of radiation transmission, a series of complicated changes take place in the light after sunlight passes through the atmosphere and the water column, then to the bottom of the water, before finally entering the satellite sensor. As shown in Figure 4, the satellite sensor receives the radiation signal L t mainly in four parts.
Appl. Sci. 2019, 9,4375 7 of 17 where is the apparent emissivity, is the relative distance between the sun and the earth, and is the average solar radiation value of the band with the unit ⁄ . The value of can be calculated according to the spectral response function of the satellite sensor. The calculation formula is where and are the spectral wavelengths of the sensor with the unit , and ( ) is the solar radiation value outside the atmosphere with the unit ( ) ⁄ . Equation (3) can be calculated from the solar spectrum response function provided by the World Radiation Center.
After the model parametric setting is completed, the 6S atmospheric correction model can be run to obtain the atmospheric correction parameters of , , and . Subsequently, the apparent reflectance data of the image are converted into surface reflectance data. The conversion formulas are as follows: is the surface reflectance, and , , and are the conversion parameters calculated by the 6S atmospheric correction model. The desired surface reflectance data were, thus, obtained.

Principle of bathymetry based on passive optical remote sensing
Sunlight has a perspective in a water column. The signal received by the optical remote sensor contains reflection information of the sea bottom to the sunlight, and this information is the physical basis of the passive optical remote sensing inversion of shallow water depth. In other words, the information reflected from the seabed to the optical remote sensor is the direct reflection of underwater terrain, that is, the information source of optical remote sensing bathymetry. The attenuation coefficient of sunlight in water determines the penetrable depth of the light in water. The smaller the attenuation coefficient of light in water is, the better its transparency in water is. This phenomenon is the basic principle of water depth inversion.
In the process of radiation transmission, a series of complicated changes take place in the light after sunlight passes through the atmosphere and the water column, then to the bottom of the water, before finally entering the satellite sensor. As shown in Figure 4, the satellite sensor receives the radiation signal mainly in four parts. = + + + , where is the atmospheric radiation, is the reflection radiance of the water surface, is the water column radiance, and is the bottom radiance. The radiative transfer equation of  where L p is the atmospheric radiation, L S is the reflection radiance of the water surface, L w is the water column radiance, and L b is the bottom radiance. The radiative transfer equation of electromagnetic waves in water is the physical basis of water depth inversion. Lyzenga [28,29] and Philpot [30] developed the water depth radiation transfer equation according to Beer's law.
where L = L t − L p − L s is the radiative brightness value after atmospheric correction and radiative brightness correction of the water column in deep water, L ∞ is the radiation brightness value at infinite depth, g is the two-way diffuse attenuation coefficient in the water body, z is the water depth value, and A d is the reflected energy value of the bottom.

Remote sensing bathymetric model
The radiation transfer equation and the basic principle of water depth inversion led to the rapid development of the bathymetric method, which mainly comprises the theoretical analytical model [31], single-band water depth inversion model [32,33], multi-band log-linear model [34][35][36], and logarithmic conversion ratio model [37]. Tian [4] improved the method on the basis of the classic Stumpf model to significantly improve water depth inversion ability. Considering the risk that the denominator of the Tian model may be equal to zero, we added the constant factor a to the bathymetric model. The a value may not increase model accuracy, but it can increase the stability of the model to a certain extent. The passive optical remote sensing bathymetric model used in this study was named ILCRM, expressed as follows: where Z is the water depth, R w (λ i ) and R w λ j are the reflectance data of the blue and green bands, respectively, a 0 and a 1 are regression coefficients, m and n are the adjustment factors, and a is the added constant factor (i.e., users can customize its value, such as setting it to 1.01). Because blue and green bands have strong transmission ability to water body, the signals received by sensors contain seabed information; thus, the two bands were used in this study. This method belongs to the classification of semi-empirical and semi-analytical models. A part of the water depth obtained from LiDAR waveform data was used as the control point to calculate the model parameters of a 0 , a 1 , m, and n. The other part of the water depth data was used as the validation point to evaluate the accuracy of the bathymetric model.

Depth extraction from waveform data
Extracting the signals of the water surface and the bottom part from an echo waveform is the key to retrieving bathymetry based on the LiDAR waveform. The retrieving methods and parameters of the sounding system were not disclosed, but they seem to apply the commonly used peak detection method [9,18] and mathematical fitting method [12,38,39]. A burr phenomenon always exists in the echo waveform signals arising from noise. A big error occurs when the peak detection method is directly used to extract the time in which the echo signals return to the water surface and the bottom. The method of filtering, denoising, and processing affects the response of the water component to the seabed waveform. Although the water component is not the focus of our extraction, it has great influence on the seabed waveform. By analyzing the radiation transmission of LiDAR in water bodies [11], we find that the received echo signals approximate a Gaussian distribution [21,39]; thus, we use the multi-Gaussian function fitting to extract the target echo signal.
Appl. Sci. 2019, 9,4375 9 of 17 where f (t) is the fitting curve expression based on the Gaussian function, a i , b i , and c i are the parameters of the function, i is the number of Gaussian function, d is the correction parameter, and t is time.
We construct the objective function as where t beg is the start time of the echo signal reception, and t end is the end time of the echo signal reception. The Levenberg-Marquardt algorithm is used to solve the optimization problem as follows: Given a reasonable initial value, we can derive the coefficient a i , b i , c i d after optimal operation and then use the peak detection method to extract the peak position of the fitting curve for the echo signals of the water surface and the bottom part as (t 1 , p 1 ) and (t 2 , p 2 ), respectively. The expression of the peak detection method is shown in Equation (12). The difference between the x-coordinates of the two peak values is the time difference of the return of echo signal from the water surface and the bottom. The bathymetry can be calculated as shown in Equation (13).
where f (t) is the sum of the Gaussian function used to iterate and fit the LiDAR echo signal, di f f (·) is the difference between two consecutive points, i.e., di f f ( f (t)) = f (t + 1) − f ; sign(·) is a sign function to depict whether the input value is negative, positive, or zero, in which the return value is −1, 1, or 0, respectively, f ind(·) returns the indices of the elements to satisfy the expressions inside it, c w is the speed of light in water, and θ w is the local incidence angle in water. The blue points in Figure 5 represent the echo signals of LiDAR, while the red curve represents the fitting curve based on the Gaussian function. where is the start time of the echo signal reception, and is the end time of the echo signal reception.
The Levenberg-Marquardt algorithm is used to solve the optimization problem as follows: = ( , , ) .
(11) Given a reasonable initial value, we can derive the coefficient , , after optimal operation and then use the peak detection method to extract the peak position of the fitting curve for the echo signals of the water surface and the bottom part as ( , ) and ( , ), respectively. The expression of the peak detection method is shown in Equation (12). The difference between the x-coordinates of the two peak values is the time difference of the return of echo signal from the water surface and the bottom. The bathymetry can be calculated as shown in Equation (13). = ( ( ( ( ( )))) 0) + 1, is a sign function to depict whether the input value is negative, positive, or zero, in which the return value is −1, 1, or 0, respectively, (•) returns the indices of the elements to satisfy the expressions inside it, is the speed of light in water, and is the local incidence angle in water. The blue points in Figure 5 represent the echo signals of LiDAR, while the red curve represents the fitting curve based on the Gaussian function.

Datum conversion
The bathymetry obtained by resolving the LiDAR waveform corresponds to instantaneous water depth. The bathymetric data need to be corrected by using tide data to derive the water depth value for the acquisition time of the optical remote sensing image. The relationship of instantaneous sea surface, mean sea level, depth, and water bottom is shown in Figure 6. However, the LiDAR data were acquired from 11:00 a.m. to 2:00 p.m. on 19 January 2013, and the collection time of each point differed. Therefore, performing a tide correction for each control point and validation point is necessary to eliminate the errors caused by different tides.

Datum conversion
The bathymetry obtained by resolving the LiDAR waveform corresponds to instantaneous water depth. The bathymetric data need to be corrected by using tide data to derive the water depth value for the acquisition time of the optical remote sensing image. The relationship of instantaneous sea surface, mean sea level, depth, and water bottom is shown in Figure 6. However, the LiDAR data were acquired from 11:00 a.m. to 2:00 p.m. on 19 January 2013, and the collection time of each point differed. Therefore, performing a tide correction for each control point and validation point is necessary to eliminate the errors caused by different tides. Figure 5. LiDAR echo signal detection based on the fitting of Gaussian function.

Datum conversion
The bathymetry obtained by resolving the LiDAR waveform corresponds to instantaneous water depth. The bathymetric data need to be corrected by using tide data to derive the water depth value for the acquisition time of the optical remote sensing image. The relationship of instantaneous sea surface, mean sea level, depth, and water bottom is shown in Figure 6. However, the LiDAR data were acquired from 11:00 a.m. to 2:00 p.m. on 19 January 2013, and the collection time of each point differed. Therefore, performing a tide correction for each control point and validation point is necessary to eliminate the errors caused by different tides.  The LiDAR bathymetry data were calculated using the same reference surface as the optical image for the satellite-derived bathymetry. MLLW data (referring to the depth data in the Figure 6) are commonly used as the reference for bathymetric mapping. Therefore, MLLW data were selected as the starting surface. The conversion formula is as follows:

Results
Bathymetric inversion around Ganquan Island was carried out based on the LiDAR waveform data and the QuickBird image. We compared the depth detected by the multi-Gaussian function and the LiDAR processing software. Moreover, the inversion capability from combining active and passive remote sensing was evaluated from the aspects of overall accuracy and specific depth range accuracy.

The Result of Radiation Calibration and Atmospheric Correction
The water body is the most important area after radiometric and atmospheric correction. Two sites at different depths of water and one site of vegetation were selected to observe the changes of the spectral curve after radiation calibration and atmospheric corrections, as shown in Figure 7. Each row of spectral curves represents a site marked in the left image. The LiDAR bathymetry data were calculated using the same reference surface as the optical image for the satellite-derived bathymetry. MLLW data (referring to the depth data in the Figure 6) are commonly used as the reference for bathymetric mapping. Therefore, MLLW data were selected as the starting surface. The conversion formula is as follows:

Results
Bathymetric inversion around Ganquan Island was carried out based on the LiDAR waveform data and the QuickBird image. We compared the depth detected by the multi-Gaussian function and the LiDAR processing software. Moreover, the inversion capability from combining active and passive remote sensing was evaluated from the aspects of overall accuracy and specific depth range accuracy.

The Result of Radiation Calibration and Atmospheric Correction
The water body is the most important area after radiometric and atmospheric correction. Two sites at different depths of water and one site of vegetation were selected to observe the changes of the spectral curve after radiation calibration and atmospheric corrections, as shown in Figure 7. Each row of spectral curves represents a site marked in the left image.

Comparison of Different Waveform Processing Means
The LiDAR waveform data were fitted to extract the water depth using the multi-Gaussian function. The bathymetric map of the whole water area (Figure 8) was obtained after tidal correction and data-level conversion.
Appl. Sci. 2019, 9,4375 11 of 17 The LiDAR waveform data were fitted to extract the water depth using the multi-Gaussian function. The bathymetric map of the whole water area (Figure 8) was obtained after tidal correction and data-level conversion. The water depth profile of northwest Ganquan Island is shown in Figure 9. The horizontal coordinate represents the offshore distance, and the vertical coordinate represents the water depth value. In the figure, the yellow points represent the sea level, and the green and blue points represent water depth at different ranges. The slope gradient was gentle, and the depth changed slightly from 3 to 17 m. Validation points were selected from the profile at intervals of 10 m to quantify the relationship between the depth inversion result of the multi-Gaussian function and the officially published value by the LiDAR processing software. The scatter plot of the two depth values is shown in Figure 10. The R 2 of the result of the multi-Gaussian function could reach 0.9 with an MRE of 5.6% in all water depth ranges. However, the MRE was relatively large in extremely shallow waters, and the time difference of the returning echo signals from the water surface and the bottom part was somewhat small and difficult to distinguish. The MAE of the whole water depth was 58.2 cm. When the water The water depth profile of northwest Ganquan Island is shown in Figure 9. The horizontal coordinate represents the offshore distance, and the vertical coordinate represents the water depth value. In the figure, the yellow points represent the sea level, and the green and blue points represent water depth at different ranges. The slope gradient was gentle, and the depth changed slightly from 3 to 17 m.
Appl. Sci. 2019, 9,4375 11 of 17 The LiDAR waveform data were fitted to extract the water depth using the multi-Gaussian function. The bathymetric map of the whole water area (Figure 8) was obtained after tidal correction and data-level conversion. The water depth profile of northwest Ganquan Island is shown in Figure 9. The horizontal coordinate represents the offshore distance, and the vertical coordinate represents the water depth value. In the figure, the yellow points represent the sea level, and the green and blue points represent water depth at different ranges. The slope gradient was gentle, and the depth changed slightly from 3 to 17 m. Validation points were selected from the profile at intervals of 10 m to quantify the relationship between the depth inversion result of the multi-Gaussian function and the officially published value by the LiDAR processing software. The scatter plot of the two depth values is shown in Figure 10. The R 2 of the result of the multi-Gaussian function could reach 0.9 with an MRE of 5.6% in all water depth ranges. However, the MRE was relatively large in extremely shallow waters, and the time difference of the returning echo signals from the water surface and the bottom part was somewhat small and difficult to distinguish. The MAE of the whole water depth was 58.2 cm. When the water Validation points were selected from the profile at intervals of 10 m to quantify the relationship between the depth inversion result of the multi-Gaussian function and the officially published value by the LiDAR processing software. The scatter plot of the two depth values is shown in Figure 10. The R 2 of the result of the multi-Gaussian function could reach 0.9 with an MRE of 5.6% in all water depth ranges. However, the MRE was relatively large in extremely shallow waters, and the time difference of the returning echo signals from the water surface and the bottom part was somewhat small and difficult to distinguish. The MAE of the whole water depth was 58.2 cm. When the water depth exceeded 15 m, the MAE became larger. This phenomenon may be due to the weak echo signal from the bottom part of the deep water that is difficult to identify.
Appl. Sci. 2019, 9,4375 12 of 17 depth exceeded 15 m, the MAE became larger. This phenomenon may be due to the weak echo signal from the bottom part of the deep water that is difficult to identify.

Overall Accuracy of Inversion Result
The data from the Aquarius LiDAR waveform were solved as survey data. A bathymetric inversion of the QuickBird image was carried out. Control and validation points were selected uniformly and randomly, and their numbers were 83 and 91, respectively. The distribution of the points is shown in Figure 11a. The regression relationship between the survey data and the inversion data of the control points was established with Equation (8). The Newton iteration method was applied to accomplish the regression of model parameters. The optimal regression parameters were , , , , and the R 2 of the fitting was 0.96. Control and validation points were then used with ILCRM to retrieve the depth value. R 2 , MRE, and MAE were used as the evaluation indicators of accuracy. The overall accuracy of the inversion result is shown in Figure 12.

Overall Accuracy of Inversion Result
The data from the Aquarius LiDAR waveform were solved as survey data. A bathymetric inversion of the QuickBird image was carried out. Control and validation points were selected uniformly and randomly, and their numbers were 83 and 91, respectively. The distribution of the points is shown in Figure 11a.
Appl. Sci. 2019, 9,4375 12 of 17 depth exceeded 15 m, the MAE became larger. This phenomenon may be due to the weak echo signal from the bottom part of the deep water that is difficult to identify.

Overall Accuracy of Inversion Result
The data from the Aquarius LiDAR waveform were solved as survey data. A bathymetric inversion of the QuickBird image was carried out. Control and validation points were selected uniformly and randomly, and their numbers were 83 and 91, respectively. The distribution of the points is shown in Figure 11a. The regression relationship between the survey data and the inversion data of the control points was established with Equation (8). The Newton iteration method was applied to accomplish the regression of model parameters. The optimal regression parameters were , , , , and the R 2 of the fitting was 0.96. Control and validation points were then used with ILCRM to retrieve the depth value. R 2 , MRE, and MAE were used as the evaluation indicators of accuracy. The overall accuracy of the inversion result is shown in Figure 12. The regression relationship between the survey data and the inversion data of the control points was established with Equation (8). The Newton iteration method was applied to accomplish the regression of model parameters. The optimal regression parameters were a 0 , a 1 , m, and n, and the R 2 of the fitting was 0.96. Control and validation points were then used with ILCRM to retrieve the depth value. R 2 , MRE, and MAE were used as the evaluation indicators of accuracy. The overall accuracy of the inversion result is shown in Figure 12.  Figures 11a and 11b show the scatter plots of the control points and the validation points, respectively. The derived depths of both points were well fitted with the survey data. However, the result of the validation points was more discrete than that of the control points, and the fitting coefficients of R 2 were 0.95 and 0.96, respectively. The water depth of the inversion was in the range of 2-22 m. The overall MAEs were 0.80 m and 0.95 m, and the overall MREs were 7.0% and 8.9%, respectively.

Accuracies of the Specific Depth Range
Water depth was divided into 10 ranges at the interval of 2 m to analyze the inversion ability of different specific depth ranges, namely, 2-4, 4-6, 6-8, 8-10, 10-12, 12-14, 14-16, 16-18, 18-20, and 20-22 m. MAE and MRE were used as the approximates of evaluation. Figure 13 shows     Figures 11a and 11b show the scatter plots of the control points and the validation points, respectively. The derived depths of both points were well fitted with the survey data. However, the result of the validation points was more discrete than that of the control points, and the fitting coefficients of R 2 were 0.95 and 0.96, respectively. The water depth of the inversion was in the range of 2-22 m. The overall MAEs were 0.80 m and 0.95 m, and the overall MREs were 7.0% and 8.9%, respectively.

Discussion
Visible light is difficult to be quantitatively analyzed due to the variability and complexity of hydrological and atmospheric factors during their propagation. The improved band ratio algorithm, which is used to retrieve depth data from passive remote sensing, belongs to the classification of semi-analytical models. In this model classification, the detection capability of water depths is limited by the depth range of control points. Given that the control points provided by the Aquarius LiDAR data are in the range of 0-20 m, the derived bathymetry can be narrowed to a relatively small range.
The extraction of echo signals from extremely shallow waters and the bottom of deep waters has certain limitations. Further research is needed to identify and extract LiDAR waveform signals. The incident angle is a key parameter of bathymetric inversion. This study assumes that the incident angle for laser launching into waters is equal to the launch angle of the LiDAR. However, the assumption is inappropriate in complex sea conditions, and it needs to be gradually improved in later experiments.

Conclusions
On the one hand, bathymetry derived by LiDAR has the advantage of high maneuverability and efficiency, but the cost of using the system is high. On the other hand, the application of passive optical remote sensing entails poor accuracy, but the cost of using it is low. On the basis of the popularity of LiDAR bathymetry for depth measurements, bathymetric inversion with passive optical image became a research topic with high value. In this study, a multi-Gaussian function was used to fit the LiDAR waveform data, and the depth data were extracted for use as the measured depth data of bathymetric inversion. Then, the improved band ratio model was applied to estimate the bathymetry of the entire water around Ganquan Island within the QuickBird image range. The main conclusions are highlighted below.
On the basis of the multi-Gaussian function to fit the LiDAR waveform data, the overall MAE of the whole water depth (3-17 m) was 5.6%, the overall MRE was 58.2 cm, and the R 2 could reach 0.97. The accuracies meet the requirements of ocean surveying and bathymetric mapping.
The bathymetric data derived from the Aquarius airborne LiDAR were divided into control points and validation points. On the basis of the improved band ratio model, a bathymetric inversion was carried out with the QuickBird image. The fitting coefficients R 2 of validation points

Discussion
Visible light is difficult to be quantitatively analyzed due to the variability and complexity of hydrological and atmospheric factors during their propagation. The improved band ratio algorithm, which is used to retrieve depth data from passive remote sensing, belongs to the classification of semi-analytical models. In this model classification, the detection capability of water depths is limited by the depth range of control points. Given that the control points provided by the Aquarius LiDAR data are in the range of 0-20 m, the derived bathymetry can be narrowed to a relatively small range.
The extraction of echo signals from extremely shallow waters and the bottom of deep waters has certain limitations. Further research is needed to identify and extract LiDAR waveform signals. The incident angle is a key parameter of bathymetric inversion. This study assumes that the incident angle for laser launching into waters is equal to the launch angle of the LiDAR. However, the assumption is inappropriate in complex sea conditions, and it needs to be gradually improved in later experiments.

Conclusions
On the one hand, bathymetry derived by LiDAR has the advantage of high maneuverability and efficiency, but the cost of using the system is high. On the other hand, the application of passive optical remote sensing entails poor accuracy, but the cost of using it is low. On the basis of the popularity of LiDAR bathymetry for depth measurements, bathymetric inversion with passive optical image became a research topic with high value. In this study, a multi-Gaussian function was used to fit the LiDAR waveform data, and the depth data were extracted for use as the measured depth data of bathymetric inversion. Then, the improved band ratio model was applied to estimate the bathymetry of the entire water around Ganquan Island within the QuickBird image range. The main conclusions are highlighted below.
On the basis of the multi-Gaussian function to fit the LiDAR waveform data, the overall MAE of the whole water depth (3-17 m) was 5.6%, the overall MRE was 58.2 cm, and the R 2 could reach 0.97. The accuracies meet the requirements of ocean surveying and bathymetric mapping.
The bathymetric data derived from the Aquarius airborne LiDAR were divided into control points and validation points. On the basis of the improved band ratio model, a bathymetric inversion was carried out with the QuickBird image. The fitting coefficients R 2 of validation points and control points were 0.95 and 0.96, respectively. In the range of 2-22 m, the overall MAEs were 0.80 m and 0.95 m, and the overall MREs were 7.0% and 8.9%, respectively. In specific depth ranges, the MAEs of all the depth ranges were within the range of 0.5-1.2 m. From the aspect of accuracy of both overall and specific depth ranges, the proposed method is an effective and inexpensive way to extract water depth information, given the combination of active and passive remote sensing. At the same time, costs can be reduced, indicating its applicability in coastal engineering in the future. Suitable depth ranges need to be determined prior extracting the depth from LiDAR waveforms and combing active and passive remote sensing. Bathymetric inversion is limited in extremely shallow waters (<2 m) and deep waters. The difference in signals from the sea surface and the bottom in extremely shallow waters is somewhat small and, thus, difficult to distinguish, and the detection accuracy is low. By contrast, in deep waters, the signal from the bottom is too weak to be detected, thus resulting in poor accuracy.