Shallow Water Measurements Using a Single Green Laser Corrected by Building a Near Water Surface Penetration Model

To reduce the size and cost of an integrated infrared (IR) and green airborne LiDAR bathymetry (ALB) system, and improve the accuracy of the green ALB system, this study proposes a method to accurately determine water surface and water bottom heights using a single green laser corrected by the near water surface penetration (NWSP) model. The factors that influence the NWSP of green laser are likewise analyzed. In addition, an NWSP modeling method is proposed to determine the relationship between NWSP and the suspended sediment concentration (SSC) of the surface layer, scanning angle of a laser beam and sensor height. The water surface and water bottom height models are deduced by considering NWSP and using only green laser based on the measurement principle of the IR laser and green laser, as well as employing the relationship between NWSP and the time delay of the surface return of the green laser. Lastly, these methods and models are applied to a practical ALB measurement. Standard deviations of 3.0, 5.3, and 1.3 cm are obtained by the NWSP, water-surface height, and water-bottom height models, respectively. Several beneficial conclusions and recommendations are drawn through the experiments and discussions.


Introduction
Airborne LiDAR bathymetry (ALB) is an accurate, cost-effective, and rapid technique for shallow water measurements [1].In general, integrated infrared (IR) and green ALB systems simultaneously and collinearly emit green and IR lasers.These two lasers are used for water bottom and water surface detection, respectively, because of their wavelength and attenuation [1][2][3].
A typical waveform of a green laser mainly includes the returns from the air-water interface, water volume backscatter, and water bottom [1,4].If the returns in the waveform are accurate, then the green laser alone may be used in an ALB system, thereby making this system compact and convenient.However, Guenther et al. (2000) rejected such a possibility due to the water surface uncertainty of green lasers [1,3].In the waveform of green lasers, the surface return is a linear superposition of the energy that is reflected from the actual air-water interface, as well as the energy backscattered from particulate materials in the water volume just under the interface [1].Therefore, the first return cannot exactly represent the water surface but reflects a certain level of penetration into the water column [5].This phenomenon is called the water surface uncertainty of the green laser [1].Such penetration is referred to as near water surface penetration (NWSP) [5].Due to the existence of the surface uncertainty, an additional IR laser is used in the integrated IR and green ALB systems to determine the actual water surface height and accurate water depth [2].The integrated IR and green (i.e., double laser) ALB systems, such as the Optech Coastal Zone Mapping and Imaging LiDAR (CZMIL), AHAB HawkEye II and HawkEye III, improve the measurement accuracy, but add to the costs and weights of ALB systems due to the double laser.To make the systems compact and cost-effective, single green laser ALB systems no longer use the primary IR laser, but only emit and receive the green laser.These systems have been extensively tested in coastal and riverine measurements, such as Fugro LADS LADS-MK3, Optech Aquarius, USGS Experimental Advanced Airborne Research LiDAR (EAARL) and EAARL-B, and RIEGL VQ-820-G [6][7][8][9][10][11][12][13][14][15].As they ignore the water surface uncertainly, the simplified systems conduct the measurements despite the loss of accuracy [5,16], and cannot comply with the requirements of highly accurate applications.
If NWSP is accurately estimated, then the surface location detected by the surface return of green laser can be corrected, and the accurate water bottom location and water depth can be determined by combining the bottom return of the green laser [1].Therefore, a green ALB system can obtain accurate measurements.At present, theoretical analysis [2,17] and statistical analysis [5,18] are used to estimate NWSP.Guenther (1985) theoretically analyzed the generation mechanism of NWSP, applied LiDAR equation to estimate the time difference between the air-water interface return and water column return, and determined the relationship between time difference and NWSP.Time difference is undoubtedly related to hydrologic conditions and ALB system parameters [1].The theoretical analysis provides the NWSP estimation and factors that influence NWSP.However, the transmission beam is assumed to be a triangular pulse characterized by multiple forward downward scattering, high-angle backscattering, and multiple forward scattering back to the surface [17].In general, the assumption is inconsistent with the actual situation, such as the Gaussian pulse used in Optech CZMIL [19].In addition, many parameters used in the LiDAR equation should be manually given.Hence, the theoretical analysis method should be improved further.Statistical analysis provides another method to estimate NWSP.Mandlburger et al. (2013) used the water surface height derived from the IR laser as a reference to statistically analyze the distributions of NWSP in ponds and rivers.The results show that NWSP ranges from 10 to 25 cm, and is high in turbid water.Moreover, it is possible to approximate the reference water surface using green surface returns by statistical distribution correction [5].The statistical analysis method is simple and efficient in calm and clear rivers and lakes, and can accurately reflect the spatial distribution of NWSP.However, this method may be affected by the undulating water surface and turbid coastal waters.
The aforementioned contributions promote our understanding of NWSP and its influencing factors.However, the precise relationship between NWSP and its influencing factors has yet to be established.Therefore, the current study proposes a method to develop the relationship model and conduct the shallow water measurements using a corrected single green laser.The structure of this paper is as follows.Section 1 introduces the existing ALB systems and studies on NWSP.Section 2 provides the detailed method of building the NWSP model.Section 3 deduces the mathematical height models of the water surface and water bottom points.Section 4 presents the validation and analysis of the proposed method through experiments.Section 5 provides the corresponding discussion.Lastly, Section 6 presents several beneficial conclusions and recommendations that are drawn from the experiments and discussions.

Comprehensive NWSP Model
Guenther et al. (1985) used the LiDAR equation to estimate the time difference between the air-water interface and the water column returns.The time difference is associated with NWSP, as well as relating to the water turbidity, transmitted pulse width, pulse detection algorithm, beam scanning angle (i.e., angle of the incident beam to the local vertical), and laser spot diameter on the water surface [1].Among these factors, the transmission pulse width of an ALB system is fixed, and the performance of the pulse-detection algorithm adopted in the waveform detection is known.The two effects can be regarded as constants.The hydrological parameters of seawater, such as turbidity; and the measurement parameters of the ALB systems, such as the beam scanning angle and spot diameter, mainly influence NWSP.Therefore, NWSP (i.e., ∆d) can be expressed as follows: where β 1 -β 7 are the model coefficients; ϕ and F are the scanning angle and spot diameter of the laser, respectively; and Turb is the turbidity of the measured sea water.
The correlations between the turbidity and suspended sediment concentration (SSC) have been investigated through extensive experiments [20][21][22][23][24][25][26].Although turbidity depends on SSC, as well as particle composition and size distribution, many experiments have shown that a good linear relationship exists between turbidity and SSC [20,22,[24][25][26].Therefore, Turb is expressed as follows: where C is SSC, and a and b are the coefficients that vary with different regions and time.However, these coefficients can be regarded as constants in the same region and in a short period [20,22,[24][25][26].
F can be calculated using the sensor height H and LiDAR divergence angle γ [27] using the following equation: The ALB system determines the transmitted pulse characteristics (i.e., initial radius and divergence angle) [28].Thus, γ should no longer be included in the model as an independent variable.We obtain the following equation when Equations ( 2) and (3) are substituted into Equation (1): Equation ( 4) is also the NWSP model proposed in this study.By employing the NWSP (i.e., ∆d 0 ), which is obtained by subtracting the known water surface height from that of the green laser, as a reference, the deviation (i.e., ε ∆d ) of the NWSP model can be calculated using the following equation: where h g s is the water surface height derived from the pulse waveforms of the green laser.h 0 s is the known water surface height, which can be obtained using an IR laser scanner mounted on the same platform as the green laser scanner; this setup is called the Mandlburger's calibration scheme [5].

Development of the NWSP Model
The NWSP model presented in Equation (4) can be developed using the following ALB measurement and hydrological data: (1) Measurement data of the ALB system, such as the water surface height derived from the green laser h g s and that derived from IR laser h r s , beam scanning angle ϕ, and sensor height H.
(2) Hydrological data of the measured water area, such as the SSC of the surface layer C.
In the preceding data, h g s , h r s , ϕ, and H can be extracted from the ALB records; and C can be obtained through the on-site sampling and laboratory analysis.After obtaining these data, the NWSP equation presented in Equation ( 4) can be established in each of the water surface points.The matrix form can be expressed as follows: V n×1 = B n×7 X 7×1 − l n×1 (6) where n is the point number.
X can be calculated based on the least squares principle as follows:

Variable Selection of the NWSP Model
Theory and experience give only a general direction as to which of a pool of candidate variables (including transformed variables) should be included in the regression model.The actual set of predictor variables used in the final regression model must be determined by analysis of the data.Determining this subset is called the variable selection problem.The goal of variable selection becomes one of parsimony: achieving a balance between simplicity (i.e., as few regressors as possible) and fit (i.e., as many regressors as needed) [29].The variable selection of the regression model can be determined using a stepwise regression procedure.
The NWSP model shown in Equation ( 4) is only an initial model that considers the main influencing factors and can be optimized by stepwise regression.Moreover, t-test is adopted to conduct significance tests on the regression coefficients of the NWSP model.Detailed descriptions of stepwise regression and t-test theories are appended in Appendix A. r air = (sin ϕ, cos ϕ), r water = (sin θ, cos θ) (8

Height Models of the Green ALB System
where O is the scanner origin, r air and r water are the unit vectors of the laser propagating in the air and water, respectively; ∆t air and ∆t water are the corresponding round trip times; and c air and c water are the corresponding velocities, respectively.Moreover, ϕ and θ are the incidence angle in the air and refracted angle in the water, respectively.If the IR and green lasers are used for the measurement, then Equation (8) can be rewritten as follows: where t 0 -t 3 are depicted in Figure 1.
Remote Sens. 2017, 9, 426 5 of 18 Once t 1 is replaced by the round trip time of the green surface return t 2 , the coordinates of the water surface and water bottom points can be achieved only by the green laser.In the substitution, it is the key to determining the accurate time delay ∆t 12 of the green surface return relative to the IR surface return.

(
) ( ) where t0-t3 are depicted in Figure 1.Once t1 is replaced by the round trip time of the green surface return t2, the coordinates of the water surface and water bottom points can be achieved only by the green laser.In the substitution, it is the key to determining the accurate time delay Δt12 of the green surface return relative to the IR surface return.
Δt12 is associated with NWSP. Figure 1 shows that Δt12 can be expressed as follows: Thereafter, we obtain the following equation: Equation (9) shows that the vector model m of the sea surface point can be calculated as follows: ∆t 12 is associated with NWSP. Figure 1 shows that ∆t 12 can be expressed as follows: Thereafter, we obtain the following equation: Equation (9) shows that the vector model m of the sea surface point can be calculated as follows: where n is the vector model of the water surface point obtained by the green surface return.
The vector model r of the water bottom point can be calculated as follows: where u is the vector model of the water bottom point and is obtained by the green bottom return, and e is its bias.By combining Equations ( 13) and ( 14), the coordinates of the water surface and water bottom points relative to the scanner origin O can be respectively presented as follows: where (S s , h s ) are the coordinates of the water surface point relative to O, and (S b , h b ) are those of the water bottom point.The other symbols are the same as those in the previous equations.
If we consider only the height, then h s and h b can be obtained by combining the green height and the estimation of the NWSP model: sin 2ϕ (17) where h g s and h g b are the heights of the water surface and water bottom points, respectively, which are derived from the green laser.∆d model is the estimation of the NWSP model.
If we have an IR laser from the IR scanner mounted on the same platform as the green laser scanner, we can use the heights derived from the IR laser as reference.The errors of h s and h b (i.e., ε s and ε b ) can be expressed as follows: where h r s and h rg b represent the heights of the water surface and water bottom points, respectively, which have been derived from the IR and green lasers.

Data Acquisition
To validate the reliability and accuracy of the proposed method, an ALB measurement was conducted in December, 2014, in a high turbidity water area (i.e., 18 km × 6 km) near Lianyungang, Jiangsu Province, China.In this measurement, the weather was sunny, the sea surface was calm, and wind speed was below 10 km/h.Accordingly, the ALB data were collected using CZMIL (see the primary technical parameters in Table 1).The CZMIL system emits green (λ = 532 nm) and IR (λ = 1064 nm) laser pulses in a collinear manner [19].The beam scanning pattern is circular with a fixed angle [30].The diffuse attenuation coefficient (Kd) derived from the moderate resolution imaging spectroradiometer (MODIS) was above 1.5 m −1 .Prior to the measurement, the ground control points were set to calibrate the CZMIL system.During the measurement, five SSC sampling stations were arranged in the measurement area (see in Figure 2).Seawater samples were collected at 0.5 m below the water surface of each sampling station in situ using horizontal water samplers and submitted to the laboratory for analysis.Each water sample was filtered, dried, and weighed.Table 2 lists each station's SSC.After the field measurement of 16 lines, the CZMIL data were processed using the Optech HydroFusion software, and the 3D point cloud data of the water surface and water bottom points were obtained.Figure 2 shows the data.

Construction and Optimization of the NWSP Model
By considering the SSC sampling station as the center, a representative water area of 100 m × 100 m is selected.A total of 16,076 water surface point pairs of IR and green lasers are obtained in the five representative water areas.By using the water surface height of the IR laser as reference, NWSP Δd in each point pair can be calculated using Equation (5).In addition, the sensor height H and beam scanning angle φ in each measuring point are extracted from the raw ALB records, and the SSC of the surface layer C in each sampling station is also obtained from Table 2. Table 3 lists the statistical results of the four-type data.Evidently, Δd ranges from 17.2 cm to 39.4 cm and has a mean of 28.6 cm and standard deviation of 3.8 cm.This change implies that the effect of NWSP on the water surface height derived from the green laser is significant.Thus, an NWSP model should be built to improve the water surface height of the green laser.

Construction and Optimization of the NWSP Model
By considering the SSC sampling station as the center, a representative water area of 100 m × 100 m is selected.A total of 16,076 water surface point pairs of IR and green lasers are obtained in the five representative water areas.By using the water surface height of the IR laser as reference, NWSP ∆d in each point pair can be calculated using Equation (5).In addition, the sensor height H and beam scanning angle ϕ in each measuring point are extracted from the raw ALB records, and the SSC of the surface layer C in each sampling station is also obtained from Table 2. Table 3 lists the statistical results of the four-type data.Evidently, ∆d ranges from 17.2 cm to 39.4 cm and has a mean of 28.6 cm and standard deviation of 3.8 cm.This change implies that the effect of NWSP on the water surface height derived from the green laser is significant.Thus, an NWSP model should be built to improve the water surface height of the green laser.The four-type data of 16,076 point pairs are used for statistical analysis to assess the change of ∆d with the beam scanning angle ϕ, sensor height H and SSC of the surface layer C. Figure 3 lists the statistical relationships.Evidently, the relationships of ∆d varying with ϕ and H are approximately linear (see Figure 3a,c), and ∆d changing with C is nonlinear (see Figure 3e).These relationships show that the proposed empirical model of ∆d (see Equation ( 4)) is reasonable.A total of 14,290 of the point pairs are used to build the NWSP model shown in Equation ( 4) by Equations ( 6) and (7).The rest are applied to test the model.The model coefficients β1-β7 are estimated via the regression analysis (see Table 4).To guarantee that these coefficients are significant in the model, t-test is adopted to conduct the hypothesis tests on these regression coefficients.The standard error (SE), t-statistic (t), and p-value (p) of each coefficient is calculated (see Table 4).Accordingly, comparing the p-values of these coefficients shows that all coefficients,  4) by Equations ( 6) and (7).The rest are applied to test the model.The model coefficients β 1 -β 7 are estimated via the regression analysis (see Table 4).To guarantee that these coefficients are significant in the model, t-test is adopted to conduct the hypothesis tests on these regression coefficients.The standard error (SE), t-statistic (t), and p-value (p) of each coefficient is calculated (see Table 4).Accordingly, comparing the p-values of these coefficients shows that all coefficients, except β 5 and β 6 , are higher than the standard α = 0.05 cutoff, thereby indicating multicollinearity in the comprehensive model shown in Equation ( 4) and that the model cannot substantially represent NWSP.To solve this problem, stepwise regression is adopted to optimize the model.Equation (19) shows the optimized NWSP model.Relative to Equation ( 4), the optimized model drops ϕ 2 and H, thereby indicating that the impact of the beam scanning angle ϕ on NWSP can be described through a simple linear function.Upon optimization, the p-values of the remaining parameters in the optimized model are below α.Therefore, the parameters in Equation ( 19) are statistically significant and should be included in the model.To analyze the contribution of each parameter in the optimized NWSP model, the standardized coefficients β of each parameter in the model are calculated (see Table 5).A detailed description of standardized coefficients theory is appended in Appendix A. The contribution of the surface layer SSC is the largest, followed by the beam scanning angle ϕ and ALB sensor height H.To accurately analyze the effects of the preceding parameters on NWSP, Equation ( 19) and these coefficients in the optimized NWSP model are used to reflect the relationships of ∆d varying with ϕ, H, and C. Figure 3 shows these relationships.Relative to the statistical relationship, the relationship between each parameter and ∆d in the optimized model is considerably clear.The impact of ϕ on ∆d is a linear positive correlation (see Figure 3b), that is, ∆d increases with ϕ.The impact of H on ∆d has a monotone decreasing trend (see Figure 3d).In a relatively small range of C, ∆d increases with C, whereas ∆d decreases gradually with C when C increases to a certain range (see Figure 3f).The phenomenon of ∆d varying with C is consistent with reality.For an ALB system with a specific laser frequency or wavelength, C directly determines the surface density of seawater, thereby affecting NWSP based on Lambert's law [31].When the value of C is high, the laser penetration ability will decrease, and ∆d will inevitably decrease.In terms of the CZMIL instruments, the turning point of NWSP changing with the SSC of the surface layer appears at C as 220 mg/L (see Figure 3f).

Height Calculation
After obtaining the NWSP model, the heights of the water surface and water bottom points in the measured water area can be calculated using only a green laser and the NWSP model.The process is as follows: First, the surface layer SSC at any location of the measured area is determined by an inverse distance weighting (IDW) interpolation using SSC and the coordinates of the SSC sampling stations.IDW is extensively used to interpolate the spatial data (e.g., SSC) and considered an effective method [32][33][34].A detailed description of IDW theory is appended in Appendix A. The interpolation adopts the assumption that the change of the surface layer SSC among these sampling stations is gradual.
Second, NWSP ∆d, is calculated by Equation ( 19) using ϕ, H, and C in each of the water surface points of the green laser.
Third, the heights of the water surface and water bottom points h s and h b are calculated using Equation (17).
Lastly, the accuracy of the height is assessed by comparing the calculation height with the reference.In the comparison, the reference heights of the water surface and bottom points h r s and h rg b are provided by the integrated IR and green lasers.

Accuracy Analysis for the NWSP Models
Using the data of the remaining 1786 point pairs in the representative water areas, the comprehensive NWSP model shown in Equation ( 4) and optimized model shown in Equation ( 19) are evaluated.The ∆d in each point pair can be obtained by considering the beam scanning angle, sensor height and surface layer SSC in the two models.The NWSP model errors can be calculated by comparing the two NWSP models with the height difference derived from the IR and green surface returns.To perform statistical analysis for the model errors, the statistical parameters and corresponding probability density function (PDF) curves of the two model errors are listed in Table 6 and drawn in Figure 4.The statistical parameters and distributions of the two model errors are nearly the same.The statistical result shows that the two models are nearly equivalent in accurately reflecting NWSP.Hence, the optimized model simplifies the expression and improves the accuracy of the comprehensive model.The preceding analysis also demonstrates that the modeling method proposed in this study is efficient.4) and optimized model shown in Equation ( 19) are evaluated.The Δd in each point pair can be obtained by considering the beam scanning angle, sensor height and surface layer SSC in the two models.The NWSP model errors can be calculated by comparing the two NWSP models with the height difference derived from the IR and green surface returns.To perform statistical analysis for the model errors, the statistical parameters and corresponding probability density function (PDF) curves of the two model errors are listed in Table 6 and drawn in Figure 4.The statistical parameters and distributions of the two model errors are nearly the same.The statistical result shows that the two models are nearly equivalent in accurately reflecting NWSP.Hence, the optimized model simplifies the expression and improves the accuracy of the comprehensive model.The preceding analysis also demonstrates that the modeling method proposed in this study is efficient.The heights of the water surface and bottom points derived from the integrated IR and green lasers are used as references.The heights obtained by the height models depicted in Equation ( 17) are used to compare the reference heights.For the statistical analysis of the height errors, Table 7 lists the statistical parameters of the model errors and Figure 5 shows the corresponding PDF curves.Evidently, the accuracy of water bottom height is higher than that of the water surface height.Approximately 92.5% of the water surface height errors are below 10 cm, whereas 82.1% of the water bottom height errors are under 2 cm.The reason is that the accuracy of the surface layer SSC possibly affects the estimation of the water surface height.In the surveyed water area, the SSC evidently changes from 110 mg/L to 315 mg/L (see Table 2).The IDW interpolation may decrease the SSC accuracy at the SSC abrupt variations.The effect of SSC on NWSP is significant; thus, an inaccurate SSC will result in an inaccurate NWSP and water surface height.Nonetheless, the method proposed in this study achieves a standard deviation of 5.3 cm in the surface height estimation and that of 1.3 cm in the water bottom height estimations, thereby proving the good performance of the proposed method in the surveyed water area.The preceding statistical results also show that the water surface height and water bottom height can be determined accurately using only the corrected green laser.estimation and that of 1.3 cm in the water bottom height estimations, thereby proving the good performance of the proposed method in the surveyed water area.The preceding statistical results also show that the water surface height and water bottom height can be determined accurately using only the corrected green laser.

Discussion
The proposed method provides a good technique to obtain accurate water surface height and water bottom height using a corrected single green laser.The following main factors determine the applications and accuracies of the proposed method.
(a) Reference water surface height Only when the actual or reference water surface height is known can the actual NWSP of the green laser be calculated, and the NWSP model and height models be established.Therefore, the reference water surface height is the basis of the proposed method.In this study's experiments, the reference height is provided by the IR laser because the double-laser ALB system (i.e., CZMIL) was adopted in the measurement.Although the processing will not influence the research result, it may

Discussion
The proposed method provides a good technique to obtain accurate water surface height and water bottom height using a corrected single green laser.The following main factors determine the applications and accuracies of the proposed method.
(a) Reference water surface height Only when the actual or reference water surface height is known can the actual NWSP of the green laser be calculated, and the NWSP model and height models be established.Therefore, the reference water surface height is the basis of the proposed method.In this study's experiments, the reference height is provided by the IR laser because the double-laser ALB system (i.e., CZMIL) was adopted in the measurement.Although the processing will not influence the research result, it may become meaningless in the actual application of the proposed method, because the aim is to use a single green laser ALB system and not a double-laser one.To address this issue, two methods are provided as follows.
(1) Addition of an assistant IR laser scanner If an assistant IR laser scanner is mounted on the same platform as the green laser ALB system, the assistant IR laser can provide the reference water surface height.The assistant laser scanner is only used prior to the measurement.This scheme was also proposed by Mandlburger et al. [5] and validated to be efficient in various water areas.
(2) Water level of calm water If the measured water area (e.g., lakes) is calm, then the reference water surface height can be determined through the water level.The water level can easily be determined, such as through the use of Global Positioning System (GPS) or leveling measurement.
(b) SSC SSC is the primary influencing factor in the NWSP model.However, the hydrological parameter can only be obtained in a few limited sampling stations.An interpolation has to be performed to calculate SSC at any location surrounded by these sampling stations.To guarantee the interpolation accuracy, the following principles need to be followed.
(1) Set the sampling stations with a certain density to ensure that the SSC from these stations can reflect the SSC variations in the water area measured.(2) Set the stations at the representative locations.Only a few stations are arranged in the water area with small SSC change.By contrast, numerous stations have evident SSC change.
(c) Wave effect As shown in Figure 6, the beam scanning angle ϕ is the angle of the incident beam to the local vertical.The angle between the incident laser and normal of the water surface at the laser point is known as the angle of incidence.If a wave slope γ exists, which is the angle of the normal to the local vertical, then the actual incident angle of the laser beam is equal to the sum of ϕ and γ.
Remote Sens. 2017, 9, 425 13 of 18 performed to calculate SSC at any location surrounded by these sampling stations.To guarantee the interpolation accuracy, the following principles need to be followed.
(1) Set the sampling stations with a certain density to ensure that the SSC from these stations can reflect the SSC variations in the water area measured.(2) Set the stations at the representative locations.Only a few stations are arranged in the water area with small SSC change.By contrast, numerous stations have evident SSC change.
(c) Wave effect As shown in Figure 6, the beam scanning angle φ is the angle of the incident beam to the local vertical.The angle between the incident laser and normal of the water surface at the laser point is known as the angle of incidence.If a wave slope γ exists, which is the angle of the normal to the local vertical, then the actual incident angle of the laser beam is equal to the sum of φ and γ.Most wave energy is typically concentrated in wind wave.Table 8 shows the wind speed and its resulting wave height H and wavelength L [35].As shown in Figure 6, the approximate range of the wave slope can be estimated using the following equation: The optimized NWSP model in Equation ( 19) is used to estimate the effect of wave on NWSP.The effect ε d can be expressed as follows: where ∆d 1 is the NWSP estimation obtained by ignoring the wave slope, whereas ∆d 2 is that obtained by considering the wave slope.
Figure 7 shows the variations of ε d with the wave slopes.Evidently, ε d ranges from −2.7 cm to 2.7 cm at wind speed of 20 km/h.When the wind speed reaches 30 km/h, ε d changes from −3.9 cm to 3.9 cm.ε d also varies with the wind speed.In our experiment, the wind speed is below 10 km/h, and the effect of the wave on the NWSP changes from −1.2 cm to 1.2 cm, which can be ignored.However, if the wind speed is above 30 km/h, then ε d changes from −4.0 cm to 4.0 cm and the maximum absolute effect is more than the accuracy of the NWSP model.Therefore, the wave effect should be considered by adding the wave slope to the beam scanning angle.

(d) Alternative methods
Shipboard echo sounding is a traditional and accurate bathymetry method.However, this method is constrained by its high operating cost, inefficiency, and inapplicability in shallow waters [36].By comparison, remote sensing methods provide substantially flexible, efficient, and cost-effective means over broad areas.The remote sensing of bathymetry is broadly categorized into two: active non-imaging and passive imaging methods [37,38].The active non-imaging method (as typified by ALB) can produce accurate bathymetric information over clear waters at a depth of up to 70 m [37].
Remote Sens. 2017, 9, 425 14 of 18 Figure 7 shows the variations of εd with the wave slopes.Evidently, εd ranges from −2.7 cm to 2.7 cm at wind speed of 20 km/h.When the wind speed reaches 30 km/h, εd changes from −3.9 cm to 3.9 cm.εd also varies with the wind speed.In our experiment, the wind speed is below 10 km/h, and the effect of the wave on the NWSP changes from −1.2 cm to 1.2 cm, which can be ignored.However, if the wind speed is above 30 km/h, then εd changes from −4.0 cm to 4.0 cm and the maximum absolute effect is more than the accuracy of the NWSP model.Therefore, the wave effect should be considered by adding the wave slope to the beam scanning angle.

(d) Alternative methods
Shipboard echo sounding is a traditional and accurate bathymetry method.However, this method is constrained by its high operating cost, inefficiency, and inapplicability in shallow waters [36].By comparison, remote sensing methods provide substantially flexible, efficient, and cost-effective means over broad areas.The remote sensing of bathymetry is broadly categorized into two: active non-imaging and passive imaging methods [37,38].The active non-imaging method (as typified by ALB) can produce accurate bathymetric information over clear waters at a depth of up to 70 m [37].Relative to ALB, the passive imaging method provides an economical and flexible method to obtain bathymetric data from space.Dörnhöfer et al. ( 2016) evaluated the potential of Sentinel-2 to retrieve water depth [39].They learned that the retrieved water depths were highly correlated with echo sounding data (r = 0.95, residual standard deviation = 0.12 m) of up to 2.5 m (Secchi disk Relative to ALB, the passive imaging method provides an economical and flexible method to obtain bathymetric data from space.Dörnhöfer et al. (2016) evaluated the potential of Sentinel-2 to retrieve water depth [39].They learned that the retrieved water depths were highly correlated with echo sounding data (r = 0.95, residual standard deviation = 0.12 m) of up to 2.5 m (Secchi disk depth: 4.2 m), even though the water depths were slightly underestimated (Root Mean Square Error (RMSE) = 0.56 m).WASI-2D is a freely available software tool for analyzing atmospherically corrected multispectral and hyperspectral imagery in both optically deep and shallow water [39,40].In deeper water, Sentinel-2A bands were incapable of allowing a WASI-2D-based separation of macrophytes and sediment, which leading to erroneous water depths [39].Relative to the expensive ALB data, the Sentinel-2A data can be obtained immediately and for free.However, the range and accuracy of the retrieved depth need to be improved further.
These methods are complementary.In deep waters, the traditional echo sounding method can be used where the remote sensing method may be inefficient.In shallow waters, the active remote sensing method can be used to perform accurate bathymetry when the traditional echo sounding is inefficient.The passive imaging method can be used when budget and accuracy demand are low.

(e) Limitations
Our experiment is conducted in a high turbidity area (Kd > 1.5 m −1 ), in which SSCs of the sampling stations range from 110 mg/L to 315 mg/L.The NWSP model is adequate to fit SSC for these SSC data.However, when SSC is extremely high, the model value may be negative; thus, the model will become inefficient.In terms of the CZMIL parameters, when SSC is above 490 mg/L, the beam scanning angle is 20 • , and the sensor height is 420 m; thus, the negative value will appear (see Figure 3f).In this case, the polynomial function in the NWSP model will become inefficient, and a Gaussian function may be appropriate to fit the extreme high SSC.

Conclusions and Suggestion
By building the NWSP model and height models, the proposed method accurately determines the water surface height and water bottom height using only a corrected green laser.Compared with the results of the green ALB system, the proposed method remarkably improves the height accuracy.Relative to the integrated IR and green ALB system, the proposed method simplifies the ALB system in size and cost.Moreover, the proposed method achieves the standard deviations of 3.0, 5.3, and 1.3 cm in the estimations of NWSP, water surface height, and water bottom height, respectively, in the experiments.
The proposed method was tested in a coastal area and further tests should be carried out in a more riverine setting.The reference water surface height is the basis of the proposed method.Therefore, an auxiliary IR laser scanner is recommended in the green laser ALB system to calibrate the system and calculate the reference water surface height.In addition, to guarantee the accuracy of the proposed method, the density and the representativeness should be considered in setting SSC sampling stations based on the measured water area.
Multicollinearity encompasses linear relationships between two or more variables [44], and can result in misleading and occasionally abnormal regression results [42].If the p-values of nearly all model parameters are greater than α, then multicollinearity may be present in the model.Thus, the model should be optimized by stepwise regression.In stepwise regression, variables are added one at a time.The order of entry of the variables is controlled by the statistics program, and the variable that will lead to the largest increase in R 2 is entered in each step.If an earlier variable becomes statistically insignificant with the addition of later variables, then such variables can be dropped from the model to prevent multicollinearity in the final regression model [42].
(3) Standardized coefficient The standardized coefficient β j of the jth variable reflects its contribution to the regression model established.Similarly, comparing the standardized coefficients of different variables can determine the priorities of their contributions to the established model and provide an easy method to determine whether one variable's effect is larger than that of another.Standardized coefficients, |β j , are beneficial when we intend to compare the relative importance of variables in the regression model [42].The larger |β j is, the more x j contributes to ∆d.The more a variable contributes to the prediction of ∆d, the more important it is [45].β j can be represented as [29,44,45]: (4) IDW method The SSC of all the sounding points can be interpolated as follows: where n is the total number of sediment sampling stations, i = 1 − n is the ith sampling station, P i is the weight of the ith sampling station, and C i is SSC of the surface layer of the ith sampling station.The weight P i is calculated as follows: where D i is the distance from the sounding point to the ith sampling station, (x, y) is the plane coordinates of the sounding point, and (x i , y i ) is the plane coordinate of the ith sampling station.

Figure 1 2 + r water ∆t water c water 2
Figure1shows the measurement mechanism of the ALB systems.The vector model for calculating the water bottom point can be expressed as follows:

Figure 1 .
Figure 1.Measurement mechanism of the airborne Light Detection And Ranging (LiDAR) bathymetry (ALB) system.The red and green colors represent the infrared (IR) and green lasers, respectively: (a) denotes the waveform detections of the IR laser and green laser; (b) shows the propagation ways of the two lasers and the bias induced by using only the green laser; (c) shows the location of the beam point P relative to the scanner origin O; t0-t3 denote the initial emission time of the laser pulse, round trip time of the IR surface return, round trip time of the green surface return and round trip time of the green bottom return, respectively; φ and θ are the incidence angle in the air and the refracted angle in the water, respectively, of the green laser; Δd is near water surface penetration (NWSP); and S and h denote the horizontal distance and vertical distance, respectively, of the beam point P relative to O.

Figure 1 .
Figure 1.Measurement mechanism of the airborne LiDAR bathymetry (ALB) system.The red and green colors represent the infrared (IR) and green lasers, respectively: (a) denotes the waveform detections of the IR laser and green laser; (b) shows the propagation ways of the two lasers and the bias induced by using only the green laser; (c) shows the location of the beam point P relative to the scanner origin O; t 0 -t 3 denote the initial emission time of the laser pulse, round trip time of the IR surface return, round trip time of the green surface return and round trip time of the green bottom return, respectively; ϕ and θ are the incidence angle in the air and the refracted angle in the water, respectively, of the green laser; ∆d is near water surface penetration (NWSP); and S and h denote the horizontal distance and vertical distance, respectively, of the beam point P relative to O.

c water 2 e
= r water ∆d cos(ϕ)c air c water − r air ∆d cos ϕ

Figure 2 .
Figure 2. Locations and scopes of the different measurements.The yellow, red, and blue colors denote the land, scope of ALB measurement, and depth contours in the measured water area.The numbers 1-5 show the locations of the five SSC sampling stations.

Figure 2 .
Figure 2. Locations and scopes of the different measurements.The yellow, red, and blue colors denote the land, scope of ALB measurement, and depth contours in the measured water area.The numbers 1-5 show the locations of the five SSC sampling stations.

1 .
Accuracy Analysis for the NWSP Models Using the data of the remaining 1786 point pairs in the representative water areas, the comprehensive NWSP model shown in Equation (

Figure 4 .Table 6 .
Figure 4. Probability density function (PDF) curves of the comprehensive and optimized near water surface penetration (NWSP) model errors.

Figure 4 .
Figure 4. Probability density function (PDF) curves of the comprehensive and optimized near water surface penetration (NWSP) model errors.

Figure 5 .
Figure 5. PDF curves of the height model errors: (a,b) show the PDF curves of the water-surface height and water-bottom height model errors, respectively.

Figure 5 .
Figure 5. PDF curves of the height model errors: (a,b) show the PDF curves of the water-surface height and water-bottom height model errors, respectively.

Figure 6 .
Figure6.Wave effect on the incident angle of the laser beam.φ is the beam scanning angle, γ is the wave slope, H is the wave height, and L is the wavelength.

Figure 7 .
Figure 7. NWSP error induced by ignoring the wave slope.The green and red dotted lines denote the error range when the wind speed is at 20 km/h and 30 km/h, respectively.

Figure 7 .
Figure 7. NWSP error induced by ignoring the wave slope.The green and red dotted lines denote the error range when the wind speed is at 20 km/h and 30 km/h, respectively.

Table 1 .
Technological parameters of the CZMIL system.

Table 2 .
Suspended sediment concentration (SSC) of the surface layer in each sampling station.Sampling Station SSC of the Surface Layer (mg/L)Remote Sens. 2017, 9, 425 8 of 18

Table 3 .
Statistical parameters of the four-type data used in the modeling.

Table 4 .
Coefficients and their significances to the comprehensive and optimized models.SE: standard error; t: t-statistics.

Table 5 .
Standardized coefficients β of the different parameters.

Table 6 .
Statistical parameters of the comprehensive and optimized near water surface penetration (NWSP) model errors.

Table 7 .
Statistical parameters of the water-surface height and water-bottom height model errors.

Table 7 .
Statistical parameters of the water-surface height and water-bottom height model errors.Statistic ParameterMax./cmMin./cm Mean/cm Std./cm

Table 8 .
Wind speed and the parameters of its resulting waves.

Speed in km/h Average Height in m Average Wavelength in m Wave Slope Range in Degrees
Wave effect on the incident angle of the laser beam.ϕ is the beam scanning angle, γ is the wave slope, H is the wave height, and L is the wavelength.

Table 8 .
Wind speed and the parameters of its resulting waves.