Kim, M.; Kopilevich, Y.; Feygels, V.; Park, J.Y., and Wozencraft, J., 2016. Modeling of airborne bathymetric lidar waveforms. *In*: Brock, J.C.; Gesch, D.B.; Parrish, C.E.; Rogers, J.N., and Wright, C.W. (eds.), *Advances in Topobathymetric Mapping, Models, and Applications*. *Journal of Coastal Research*, Special Issue, No. 76, pp. 18–30. Coconut Creek (Florida), ISSN 0749-0208.

Modeling the optical power of the lidar return waveform is performed for the radiometrically calibrated CZMIL (Coastal Zone Mapping and Imaging Lidar, Optech, Inc.) data. For this purpose, a lidar waveform simulator was developed. The theory is described based on the receiver sensitivity function, the radiative transfer equation (RTE) via the Greens function, the optical reciprocity theorem, and the small angle approximation (SAA). The SAA-based RTE is solved for the radiance distribution using the Fourier transform method. Along with the numerical algorithms, the contribution was made on the air-water and water-bottom interface peaks in the bathymetric lidar waveforms. Lacking ground truth data, a simulated waveform that best fits CZMIL data was used to estimate the optimized environmental parameters. The estimated parameters were well within the plausible natural optical properties. Compared to other approaches based on the relative intensity waveform, the simulation was applied to the absolute calibrated power. Thus, the model can be used to predict the general performance of any bathymetric lidar. This research will help design an optimized system to achieve the maximum performance. The forward modeling capability will also provide opportunities to develop advanced waveform processing algorithms, such as surface peak modeling and scattering correction. Thus, the improved quality of bathymetric lidar data contributed by this research will promote the various coastal science applications in terms of improved data accuracy and extended coverage.

## INTRODUCTION

Lidar technology has dramatically improved in the last several decades and is a main bathymetry and coastal mapping tool (Wozencraft and Millar, 2005). The primary purpose of a bathymetric lidar is to measure water depth, but it can also be used to estimate optical properties of the water, to calculate the reflectance of the sea floor, and to detect submerged objects. Thus, it is of great interest to design a bathymetric lidar system that is optimized for these multiple purposes. The system optimization requires a lidar waveform simulator. A previous effort to simulate waveforms was based on the approximation of waveforms at the phenomenal level (Mathur *et al*., 2010; Ramnath *et al*., 2010). For instance, the surface and the bottom peak was modeled as a simple Gaussian peak due to the geometric stretching of the pulse. The backscattering region was modeled as an exponential function to describe the observed shape of the waveform from the water body. The rising edge of volume backscattering near the surface and the falling edge near the bottom was modeled using error function.

While the model was a semi-empirical description of the phenomenal observation, the current effort is based on a full radiative transfer solution that emphasizes the handling of scattering dispersion. A mathematical model based on advanced radiative transfer theory can predict the waveforms upon the variations in the geometrical and optical system parameters and environmental properties. The simulation of this process is a fundamental point for the quantitative interpretation of lidar waveforms and may be described by a mathematical model (Kopilevich *et al*., 2005). The process could be modeled as shown in Figure 1.

Consider an infinitesimally thin disk of laser pulse (delta pulse) with unit radiant energy, which propagates through the coastal environment. The history of the light interaction is returned to the photo-detector as a form of optical waveform, *S*_{δ}(*t*), which is a combined response by a lidar system and the environment due to the delta pulse (Guenther, 1985). This is the most important step in the entire waveform simulation. The rest of the simulation concerns how detector-electronics respond to *S*_{δ}(*t*) . The main topic of this research is to model *S*_{δ}(*t*) using the radiative transfer solution. Since the laser pulse in the real world has finite thickness with a normalized pulse shape function, *p(t)*, the optical signal *S(t)* in radiant power can be modeled as a convolution, *S*_{δ}(*t*) * *p*(*t*).

The optical energy creates photo-electrons based on the quantum efficiency (QE) of a photo-detector. Subsequently, the photo-electrons are amplified (*e.g.*, PMT or APD) to a high current. This nonlinear amplification characteristic, α, consists of the major part of the radiometric calibration of a lidar system (Figure 10). The conversion of optical energy to electrical energy and its amplification always accompanies the time delay/stretching, which is characterized by a normalized system response function, *ω(t)*. It includes the photo-detector response (rise/transition/fall) and finite electronic system response. Thus, the photo-electronic process that creates amplified current requires another convolution as shown in Figure 1. The electric current *I(t)* is converted to voltage *V(t)* when it passes through a terminating resistor with resistance, ρ, and the voltage is digitized to digital count *D(t)* using a conversion factor, ξ [V^{−1}]. Note that the separation of the laser pulse shape *p(t)* and the function *ω(t)* is not required in this process. The total system response function as a convolution *R*(*t*) = *p*(*t*) * *ω*(*t*) is sufficient.

Since the optical signal *S*_{δ}(*t*) is the essence of the whole process, it will be described in detail. Using the above framework, a lidar waveform simulator was developed and it was validated using CZMIL data by JALBTCX. This paper consists of three sections:

(1) When a delta-pulse of a localized source propagates into an absorbing and scattering hydro-optical medium, it is constantly backscattered by the medium. It is also reflected when it encounters interfaces such as the water surface and sea floor. The theoretical basis of these processes is addressed and the general governing equation of the lidar waveform is derived.

(2) The procedures to calculate numerical solutions based on this theory is described. Since the essence of the numerical solution is the irradiance profile, detailed description is made for the properties of the irradiance function in relation to geometric and water properties to derive an efficient numerical calculation strategy.

(3) The validation of the waveform simulator is performed on the radiometrically calibrated CZMIL waveforms. The CZMIL waveforms collected from an airport runway is used first because the quality of the simulator can be best evaluated using data from a controlled target. It is followed by many waveforms collected from the shallow to the deep coastal ocean water for the validation.

## METHODS

The governing equation for the bathymetric lidar waveform is derived based on the radiative transfer of the laser pulse, assuming the propagation within a small angle is due to highly forward scattering by ocean particles. The details of the numerical calculation are presented.

### Derivation of the General Lidar Equation

Consider a lidar receiver at *z*_{0}-plane that looks downward to collect the upward return radiance that is a result of the volume backscattering of a mono-static laser beam (Figure 2). The return power *S*_{δ}(*t*) due to a delta-pulse can be modeled (Kopilevich and Surkov, 2008) as the integral of the upward return radiance weighted by the receiver sensitivity function:

**r**is the position vector on a plane (

*z*

_{0}-plane) perpendicular to the beam direction,

**n**is a unit direction vector, and Ξ

_{u}represents the upper hemispherical integral domain for a solid angle centered to

**n**. The

*I*′

*R*(0,

**r**,

**n**) is the normalized receiver sensitivity distribution function.

The upward lidar return radiance, *I* ^{up} (*t*,0,**r**,**n**), is the solution of the radiative transfer equation (RTE) for backscattered radiation, and it can be expressed with backscattered, upward radiance distribution, *I* ^{up}(*z*,**r**,**n**), at *z* via the Greens function:

**n**′ ) point source at the position,

**r**′, on the

*z*-plane to the radiance at the position

**r**on the

*z*

_{0}-plane to the direction

**n**. The upward radiance

*I*

^{up}(

*z*,

**r**′,

**n**′) is obtained using “multiple forward, single backscattering” approximation (Dolin and Levin, 1991; Monin, 1983): where

*c*is the speed of light,

*n*is the refractive index, and

*β*is the volume scattering function (VSF). The radiance distribution of the laser pulse

*I*′

_{E}(0,

**r**,

**n**) is shown as

*I*′

_{E}(0) in Figure 2.

*I*′

_{E}(

*z*,

**r**′,

**n**″) is the radiance solution on the

**-plane of the stationary RTE for downwelling radiation due to the laser source**

*Z**I*′

_{E}(0,

**r**,

**n**) . Using Equations 2 and 3 and rearranging the integral term, Equation 1 can be rewritten as

An optical reciprocity theorem is applied (Dolin and Levin, 1991; Mobley, 1994; Monin, 1983; Walker, 1994):

The square bracket in Equation 4 then becomes the “fictitious receiver radiance” *I*′_{R} (*z*,**r**′, **n**′).

A “true laser beam,” *I*_{E} (0,**r**,**n**), propagates to the *z*-plane and creates a radiance field, *I*′_{E} (*z*,**r**′,**n**″), determined by the beam diameter, divergence, and optical properties of the medium. The *I*′_{R} (*z*,**r**′, −**n**′) can be interpreted in the same manner.

A “fictitious receiver beam,” *I*′_{R} (0,**r**, −**n**), propagates to the *z*-plane and creates a radiance field, *I*′_{R} (*z*,**r** *I*′_{R}, −**n**′), determined by the receiver diameter, FOV, and medium. The “receiver beam” *I*′_{R} (0,**r**,−**n**) is, in fact, the lidar sensitivity distribution (Figure 2) to the opposite direction. Using *I*′_{R} (0,**r**,−**n**), the backscattered power of Equation 4 becomes

*E*(

**r**) ≈ ∫

*dω*(

**n**)

*I*(

**r**,

**n**) is approximated based on the small angle approximation (SAA). Equation 5 is a lidar equation for volume backscattering. For any other physical interaction, such as surface, bottom, and object return, the only change needed is to replace

*β*(

*z*,

**r**) with a proper interaction function so that Equation 5 can be used as a general lidar equation.

### Radiance Solution of SAA-based RTE

A lidar waveform can be calculated using Equation 5 with known irradiances *E*_{E,R} . The irradiance is calculated by the integral of the radiance distribution. Subsequently, the radiance distribution is obtained as a solution of the RTE:

*c*is the beam attenuation coefficient and

**ρ**= (

**r**,

*z*) is a 3D position vector. However, this 3D RTE Equation 6 does not yield an analytical solution. In order to transform 3D RTE to an approximated RTE that yields a solution, the SAA,

**n**= (

**n**

_{t,}

*n*

_{z}),

*n*

_{z}≈ 1, |

**n**

_{t}|<<1, is used based on the fact that the VSF of the natural water is strongly forward-oriented. Also, instead of regular VSF, a small angle phase function

*x*

_{s}(

*θ)*is defined within forward angles (0,

*θ*

_{s}).

The small angle scattering coefficient is defined by using the scattering coefficient *b*. Using these approximations, an SAA-based RTE is formulated as

The beam propagates along the *z*-axis. Due to the SAA, the displacement |**n**_{t} – **n**′_{t} is equivalent to the scattering angle, and the solid angle *dω*(**n**) is reduced to the infinitesimal area *d***n**_{t}. The SAA-based RTE yields an analytical solution (Walker, 1994; Bremmer, 1964; Dolin, 1964; Ishimaru, 1978). A Fourier transform of Equation 7 is performed with respect to **r** via the angular frequency conjugate **k**:

*Î*(

*z*,

**k**,

**n**

_{t}) = (2

*π*)

^{−1}∫exp[

*i*

**k**·

**r**]

*I*(

*z*,

**r**,

**n**

_{t})

*d*

**r**.

Introducing *Î*′ to Equation 8 yields an RTE of *Î*′

Taking the Fourier transform of Equation 10 w.r.t. **n**_{t} yields

The solution to the ODE Equation 11 is given by

Taking the Fourier transform of Equation 9 gives

Equating Equations 12 and 13 yields the analytical solution of the SAA-based RTE

### Irradiance Solution and the Source Functions

The irradiance is obtained via an inverse Fourier transform of the radiance from Equation 14:

Rearranging Equation 15 and using the analytical Fourier transform of 2D delta distribution yields

Plugging Equation 14 with **p** = 0 into Equation 16 yields the final irradiance equation:

Equation 17 can be rewritten as follows:

where*a*

_{s}

*= c - b*

_{s}, and

*a*

_{bs}

*(kz)*is the contribution from the small-angle scattering to attenuation of spatial harmonics of the frequency

*k*in initial distribution of a light beam radiance when the beam is propagated to

*z*(Kopilevich

*et al*., 2010; Kopilevich

*et al*., 2011).

### Calculation of the Normalized Irradiance Function

In order to describe the oblique entrance of the laser pulse and its propagation, *z – z*_{0} from Equation 18 is replaced by in-water propagation distance *h*,

*h′*(

**r**) is the in-water path length to a position

**r =**(

*r*

_{=},

*r*

_{⊥}) on the

*h*-plane. The

*r*

_{=}-axis is defined as an intersection between the

*h*-plane and a plane that includes both nadir and propagation direction. In case of an incident angle

*θ*

_{a}and corresponding refractive angle

*θ*

_{w}, the in-water distance

*h′*(

**r**) =

*h–r*

_{=}· tan

*θ*

_{w}is the function of

*r*

_{=}only. The positive

*r*

_{=}is to the direction of propagation. The

*r*

_{⊥}-axis is determined from the

*r*

_{=}-axis in the right-hand rule.

The next step is to define a Fourier transform of the source *Ĩ* in Equation 19. Consider a lidar located at altitude H. The source radiance distribution is characterized by the angle Θ (beam divergence for beam and FOV for receiver). If the in-water propagation is projected to the air, an equivalent problem is formulated. A Gaussian model and its Fourier transform are given by

*r*

_{0}(=

*H*

_{s}Θ / 2) is the radius at

*e*

^{−1}level, and for a Step-function model where

*r*

_{0}is the cutoff radius,

*J*

_{1}is the first order Bessel function of the first kind. Plugging the Fourier transform of the Gaussian source function in Equation 20 into Equation 19 yields the irradiance equation with the 0-th order first kind Bessel

*J*

_{0}where the effective in-water radius is defined as P = (Θ / 2) · (

*H*

_{s}+

*h*/

*n*) . The normalization condition of the irradiance using the angular symmetry is defined as

As the beam propagates, the forward scattering spatially redistributes radiance over the beam front. Thus, despite the dispersion due to the forward scattering, the total energy is conserved. Since the attenuation term including 1/2π is already separated, the spatial integral of the *g*_{Θ}(*h*,*r*) must be unity, which means it is the normalized irradiance function.

The normalization condition can be used to check the accuracy of the numerical calculation of *g*_{Θ} (*h*,*r*) . In order to calculate the effective attenuation *a*_{bs}*(hk)* in Equation 22, an analytical forward scattering phase function model suggested by Dolin (Dolin *et al*., 1988) is used.

Due to the oscillatory nature of the Bessel function, it is numerically integrated for the individual lobes of the Bessel function in Equation 22 with *k _{i}= J*

_{0,i}/

*r*as boundaries, where

*J*

_{0,i}is the i-th Bessel

*J*

_{0}root

*w*

_{n}for each node

### Characteristics of the Normalized Irradiance Function

The normalized irradiance function *g*_{Θ} (*h*,*r*) is the essence of waveform simulation. Therefore, it is worth investigating the characteristics of the integrand function *f(k;h,r)* in Equation 22. The set of variables used for demonstration is H = 400m, *θ*_{a} = 20^{°}, α = 7, Θ = 10mrad, and *b*_{s} = 0.3m^{−1}.

The exponential term in Equation 22 acts as a damping function that restricts the spatial frequency range of the Bessel function. The effect of h (0–20 m variation) on the two terms in the exponent is illustrated in Figure 3. Regarding the first term, *- (k*^{2}*/4)*P^{2}, the effect of *h* is virtually negligible, *h/H*_{s} ≈ 1, and P(*h*) is virtually constant. Thus, the 10 curves from varying values for *h* are tightly packed. The second term, −*a*_{bs}*(hk)h*, is the optical thickness due to the attenuation by the small angle forward scattering. Since the optical depth is proportional to geometrical distance, *h* has a significant effect, resulting in wide vertical spread from the top curve (0 m) down to the bottom curve (20 m). An additional complication is that *h* acts as a contracting factor of the *a*_{bs}*(hk)*.

This means that *a*_{bs}*(hk)* function is contracted along the spatial frequency axis for longer values of *h* and stretched for shorter values of *h*. The whole exponential term (the top group of curves in Figure 3) is the damping function in Equation 24. In the numerical calculation of the irradiance function, it determines the number of Bessel lobes to be integrated in Equation 25. The upper limit of spatial frequency *k*_{max} for the integral can be estimated from a big enough number (*e.g.*, 10).

Once the overall damping function is calculated, the most prominent feature in calculating *f(k;h,r)* is the distance *r* from the beam center because it is a contraction factor of the Bessel function. Therefore, the *f(k;h,r)* function with a bigger *r* includes a larger number of Bessel function lobes within the same spatial frequency range (Figure 4). Regarding Θ (divergence or FOV), the smaller angle (*e.g.*, Θ=2mrad) results in smaller P. Thus, the slow damping spans to relatively larger spatial frequency (Figure 4a). On the other hand, a larger angle (*e.g.*, Θ=10mrad) limits the function within small spatial frequency range (Figure 4b).

The integral of *f(k;h,r)* results in the normalized irradiance function *g*_{Θ} (*h*,*r*) . Figure 5a shows the *f(k;h,r)* for varying distances *r* from the center up to 3 m for the case of Θ=2mrad. Figure 5b shows the profile of the *g*_{Θ}(*h*,*r*) function as a result of the integral, where the curves in Figure 5a are marked as dots in Figure 5b. Also shown is *G*(*r*) = ∫*g*_{Θ}(*h*,*r*)*rdr* converging toward unity, which demonstrates that *g*_{Θ}(*h*,*r*) is accurate based on the normalization condition from Equation 23.

Regarding the scattering effect, higher scattering results in faster damping of the trailing Bessel lobes magnitude (Figure 5c).

A collection of *g*_{Θ}(*h*,*r*) profiles gives the whole picture of the dispersion due to the forward scattering. The peak of the irradiance quickly decreases with the increasing dispersion due to the forward scattering (Figure 6a). When the *g*_{Θ}(*h*,*r*) profiles are normalized by their peak values, the effect of the scattering dispersion is pronounced (Figure 6b). The irradiance profiles for Gaussian (Figure 6a) and the Step function (Figure 6c) are illustrated.

In visualizing irradiance profiles, it is useful to define an effective radius and to simplify further using the normalization condition of *g*_{Θ}(*h*,*r*) as follows:

Another intuitive, flexible definition of the radius is the distance with a certain cumulative value from the beam center. For instance, a distance with 70% of energy, *r*_{0.7}, is defined as

Using this definition, Figure 7 illustrates the irradiance dispersion pattern. The wide, white dotted lines are for the 40 mrad FOV receiver, and the narrow but fast dispersing green dotted lines are for the beam with 5 mrad divergence. Thick white and green lines represent the *r*_{0.7} defined in Equation 29 for the FOV and the beam, respectively. Thin white and green outer boundary lines represent *r*_{0.98}. Additionally, the two light blue lines defined in Equation 28 are shown so that it can be compared to *r*_{0.7}. Since the two values are defined differently, the dispersion shows a little discrepancy. However, they both serve well as an indicator of the scattering dispersion. One important thing to note is that irradiance distribution with a smaller Θ (in this case 5 mrad) disperses much quicker. The thick cyan and yellow lines represents the surface and the bottom, respectively.

### Reflection from the Surface and the Bottom

Equation 5 can be generalized for the lidar waveform for the surface and bottom reflection. Due to the non-nadir incidence to the water surface, a finite time is needed for a delta pulse to submerge. The time duration, *δt* = (2*r*_{eff} / *c*)tan*θ*_{a}, is determined by the incident angle and the effective size of the beam, *r*_{eff =} *H*_{s} Θ / 2 . The light interaction with the surface at one moment occurs only along the line segment where the delta-pulse intersects with the water surface. To derive the power from surface Fresnel reflection due to a delta pulse, , the energy during the small time duration *δt* can be equated by the 2D integral.

For any given *r*_{=} value, the practical integration range on the *r*_{⊥} -axis is finite. The increment can defined by *δr*_{=} = *cδt*/(2tan*θ*_{a}) . Cancelling *δt* from Equation 30 yields

The effective Fresnel reflectance *ρ*_{Fr} can be thought of by assuming that the distribution of the countless tiny surface wave facets acts like a Lambertian reflecting surface. In the same manner, the bottom reflection is

where the incident angle to the bottom is defined by *ϑ* = cos^{−1}(−**n** · **n _{bottom}**) with the bottom normal vector

**n**

_{bottom}and the direction of the beam propagation

**n**. In the simplest case, when the bottom surface is normal to the beam propagation plane and its slope angle is

*θ*

_{b}, then the incident angle is

*ϑ*= |

*θ*

_{w}–

*θ*

_{b}.

## RESULTS

The simulation results using the waveform simulator are compared for the CZMIL waveforms. Some of the topographic waveforms and the bathymetric waveforms from various depths are modeled by estimating the optimal environmental variables.

### CZMIL System Response Function

The CZMIL system specification is shown in Table 1. The laser pulse shape of the CZMIL laser (Bright Solution SRL) is Gaussian with 2.2 ns FWHM. CZMIL records the waveform of the laser pulse using the pick-off fiber optic cable directly hooked up to a dedicated PMT-digitizer unit. Thus, the laser pulse experiences no additional geometrical stretch, but it is only convolved with the pure electronic system response function and the resultant waveform has 2.9 ns FWHM (Figure 8a). The dots represent CZMIL data with a native 1 ns interval. Using these two observations (2.2 ns pulse duration and 2.9 ns CZMIL waveform), the pure electronic system response is calculated as a Gaussian function with 1.9 ns FWHM, which demonstrates that CZMIL has a superb detector-digitizer.

## Table 1.

CZMIL deep channel system specification.

Figure 8a is a combined result of (1) the laser pulse function, (2) the PMT characteristics, and (3) the digitizer and the rest of the electronic system. Thus, it is effectively the total system response function, which will be used as a convolution function for the delta-pulse waveform *S _{δ}*(

*t*) .

### Variability of Waveform from Different Source Functions

While typical receiver optics have sensitivity distribution characteristics of the Step function, the cross-sectional power distribution of the laser beam could be Gaussian or Step function. Consider a delta-pulse hitting the surface with a 20-degree incident angle. Due to the different intensity distribution along the cross-section, the delta-pulse return waveforms from two source functions show substantial differences (Figure 8b). The intensity of the surface return peak is related to the length of the line-segment where the surface and the slant-incident circular beam plane intersect. In case of the Step function, the line-segment is proportional to , where *r*_{0} is the effective radius of the beam. Figure 8b shows the simulated waveform exactly as predicted. However, the differences between the Gaussian and the Step function practically disappear after the overall system response function is applied to the delta-pulse return.

### QA of Simulated Return Peak from a Flat Surface

The beam divergence at *e*^{−2} level (7 mrad, Table 1) can be converted to 4.12 mrad FWHM divergence after is multiplied. When the pulse hits the ground, the temporal FWHM of the delta pulse return peak is estimated by

The predicted FWHM (4.2 ns) using Equation 33 compares well with the simulated surface return waveform from the delta-pulse (4.1 ns, Gauss_delta curve in Figure 8b). It is estimated that the 4.1 ns FWHM will be increased to 5.02 ns after the convolution with system response function and the FWHM of the simulated waveform is 5.1 ns (Gauss curve in Fig 8b). When a mono-static lidar system has FOV greater than its beam divergence, the energy of the return peak from the Lambertian surface is calculated as follows.

Using the surface reflectance ρ=0.15, Equation 34 predicts 9.21 pico-Joules upon arriving at the PMT photo-cathode out of the 3 mJ laser pulse energy. The total integrated energy from the surface return peak (Gauss curve in Fig 8b) is calculated to 9.19 pico-Joules.

The simulated waveforms are compared with CZMIL waveforms from a flat runway of the Stennis airport in Mississippi, USA. The circled area of the runway camera image in Figure 9 shows a substantial range of peak magnitudes (Figure 9a) due to the surface reflectance variation. The waveforms were normalized and the average FWHM was measured as 4.2 ns. In the previous section, the FWHM was predicted as 5.1 ns using the standard CZMIL geometry. For these runway waveforms, the sensor altitude was 387 m and the off-nadir angle was 15.4 degrees due to the platform's movement. As a result, the FWHM of the simulated waveform using proper parameters is calculated to 4.22 ns, which is a superb agreement with the real CZMIL waveform. The estimated reflectance, however, did not match well with the ground truth data.

Out of 27 waveforms, the highest and lowest two CZMIL waveforms (Hi and Lo in Figure 9a) were selected and the reflectance from the best-fit simulated waveforms were estimated. The optimized reflectance was 2.1 and 2.8%, respectively, while the average ground reflectance obtained using an ASD Fieldspec spectroradiometer (ASD Inc.) was 3.3%. The simulated waveform is calculated in terms of power. Thus, for the validation of the simulator, the CZMIL digital counts also should be converted to power (Figure 10) by applying radiometric calibration characteristics.

### Optimization of Bathymetric Waveforms

CZMIL bathymetric lidar waveforms were collected at Maui, Hawaii, on Oct 20, 2013, as a part of the National Coastal Mapping Program (NCMP) Campaign by the U.S. Army Corps of Engineers (USACE). The data collection was performed by the Joint Airborne Lidar Bathymetry Technical Center of Expertise (JALBTCX). The location of the dataset is 20.908 E, 156.698 W. The arrow in the image (Figure 11) illustrates the direction from shallow to deep water, and the figure compares 5 waveforms from varying depths (approximately 2, 10, 20, 30, and 40 m) sampled along the path. Two important features of the waveforms are found. First, the surface return peaks are co-aligned, although there are large variations in the magnitude of the surface return peaks. The large variation of surface return peaks is due to the highly variable nature of the water surface condition. Second, the overall trend of decreasing bottom return peaks demonstrates the strong correlation between the bottom return peak and the optical depth, although the bottom reflectance variation can affect the trend.

In Figure 11, each of the CZMIL waveforms from 5 different depths are over-plotted with the matching simulated waveforms. Major inversion parameters from the simulated waveform that makes the best fit to the real CZMIL waveforms are shown in Table 2. Sensor altitude and the air-incident angle variation are solely due to the aircraft kinematics. The large variation of the effective reflectance of the water surface, *ρ*_{Fr}, demonstrates the nature of complex water surface reflectance. The average volume backscattering, *β*_{π}, is higher for shallow water and lower for deep water, which is plausible because shallow water is likely to have more scattering particles. A similar trend is observed for *a+b*_{b} as well. The bottom reflectance, *ρ*_{Fr}, is independent of the depth or water optical properties.

## Table 2.

Inversion parameters from best-matching simulated waveform.

A CZMIL waveform from very shallow (2 m) water is fitted with a waveform simulator (Figure 12). The mismatch immediately after the surface peak is due to the complex vertical distribution of backscattering particles near the surface. The extended, raised tail of the bottom peak is quite different from the ideal model. When the laser beam is reflected at the bottom, the main Gaussian peak is created by the bottom-reflected radiance only within the tiny solid angle that directly goes back to the receiver, which is only 10-8 of the bottom incident energy. The majority of the energy to all other direction that was reflected during 25–35 ns must experience a short optical delay (roughly 10 ns) and then only a very tiny portion goes to the receiver.

The most likely reason for this optical delay is the high concentration of suspended mineral particles near the bottom due to the wave action in the shallow water. Another possibility could be the multiple reflections within sea grass bed. Sufficiently deep water (the 20-m waveform in Figure 12) apparently does not have the extended bottom tail.

In the case of the waveform sampled around the 11 m region, a very wide bottom peak is observed and the simulator fails to properly model with a flat bottom assumption (Figure 13). The simulator predicts a 9.2 ns FWHM for the bottom peak, but the CZMIL waveform has a much wider peak (13.8 ns). The reason for this mismatch is highly likely to be the non-flat sea floor. The bottom FWHM plot of waveforms from 5 different depths (small plot in Figure 13) demonstrates that this hypothesis is likely to be correct. The bottom peak width increases in deeper water because of geometrical stretch and beam dispersion. The only exception is the FWHM from 11 m. The average bottom plane must be slanted so that the sweeping time is significantly increased. Using 4 neighbor points, the normal vector of the bottom surface is calculated as illustrated in Figure 13 using the equation **n**_{bottom} = (**n**_{1}× **n**_{2} + **n**_{3} × **n**_{4} ) / 2, which is the average of two normal vectors (cross product). The direction of the beam **n**_{beam} is obtained via geocorrection of the lidar point cloud. Using two normal vectors, the incident angle to the sea floor are calculated and the bottom slope angle is estimated as cos^{−1}(**n**_{beam} · **n**_{bottom} ) – *θ*_{w} ≈ 14.

After applying a 14-degree slope, the predicted waveform (circle in Figure 13) matches better to the CZMIL waveform.

Simulated waveforms from 30 and 40 m are compared with CZMIL waveforms in Figure 14, where a line plot is used instead of a dot plot due to the large number of points. Deep water waveforms have a very long volume backscattering region between the surface and the bottom peaks, which gives us a better chance to extract the water's optical properties. In order to investigate the attenuation characteristics properly, a log scale plot of the 40-m waveform is illustrated (Figure 15).

When a laser pulse propagates downward, the attenuation (due to absorption and backward scattering) reduces the beam energy. Only a tiny fraction (related to *β*_{π}) of the scattered energy goes back to the receiver. The backscattered power experiences further attenuation during the return trip. This process follows the fundamental Beer's law. This simple description, however, is only a part of the total system attenuation. The combination of the FOV and the VSF creates additional attenuation. As the scattering disperses the beam, backscattering that occurs at the edge of the beam boundary is very difficult to detect because it is out of the viewing field. If the FOV is very large or scattering is vanishingly small, this additional attenuation would not be present. The waveform simulation model takes care of this effect as explained in the theory section. Thus, the biggest interest is to evaluate how a CZMIL waveform from deep water reveals its attenuating characteristics compared to the simulated waveform.

The combined total attenuation is called the “total system attenuation coefficient” (*K*_{sys}) and the factor that creates additional attenuation is called the “F-function” (Kopilevich *et al.*, 2010; Kopilevich *et al*., 2011). The red line in the main plot (Figure 15) is shown as an ideal reference that follows Beer's law, which means F-function is unity at all depths. The F-functions are calculated within a 50–300 ns region for both CZMIL data and the simulated data. The good match between the two F-functions (upper small plot in Figure 15) indicates the validity of the simulator and demonstrates the high fidelity of CZMIL's optical, PMT, and digitizer system. The lower small plot shows varying system attenuation as a function of time. In ideal modeling, the system attenuation is close to (*a+b*_{b}) and is increased with the depth due to increasing beam dispersion by forward scattering. Although it is possible to calculate a continuous *K*_{sys} profile, in practice it is highly susceptible to noise. Thus, the raw waveform was smoothed before applying a 5-m-wide moving log-linear kernel to calculate the *K*_{sys}. Although the data was smoothed first, the result (lower small plot in Figure 15) still shows large variation. In fact the high susceptibility to noise is typical in the very clear water. In turbid water *K*_{sys}, is much more stable. Despite the large variation, the result shows an overall similar pattern of gradual increase.

In typical data processing, an average *K*_{sys} is calculated within a stable volume backscattering region between the surface and bottom peaks. The waveform from 70 to 270 ns was used to calculate the average *K*_{sys} as 0.071 m-1 and the optical thickness as 1.568. The optical depths were calculated using two continuous *K*_{sys} profiles (lower small plot in Figure 15), one from the simulated waveform and the other from the CZMIL waveform. The results (1.574 and 1.529) are close to each other and are even very similar to the optical depth calculated using just two points. This demonstrates that even the simple 2-point average *K*_{sys} algorithm is good enough.

## CONCLUSIONS

A numerical calculation of the bathymetric lidar waveform is described along with the basis theory with substantial details. A pseudo-validation of the simulator was performed using radiometrically calibrated CZMIL lidar waveforms. While a ground truth for topographic data was available, the waveforms from the coastal ocean were not backed up with in situ data. Collecting a full set of optical closure data (bottom reflectance, absorption, and VSF) is difficult. Despite this limitation, it was possible to optimize the simulated waveforms to achieve highly respectable agreement with the CZMIL waveform. The water IOP's inversion results show a consistent range of values found in the typical Hawaiian water. Considering that all of the validation was performed with the waveform in optical power instead of an arbitrarily scaled waveform, the waveform simulator demonstrated a high level of calculation accuracy.

## ACKNOWLEDGMENTS

Office of Naval Research (ONR) (N00014-10-C-0085). Work by the lead author performed under USGS contract G15PC00012.

## LITERATURE CITED

*Izv. Vyssh. Uchebn. Zaved. Radiofiz*, 7, 471– 478. Google Scholar

*Izvestiya, Atmospheric and Ocean Physics*, 21, 893– 896. Google Scholar

*Reference book on the underwater vision theory*. Leningrad: Gidrometeoizdat Press, in Russian. Google Scholar

*Airborne Laser Hydrography, System Design and Performance Factors*. Rockville, Maryland: National Ocean Service, NOAA, 385 p. Google Scholar

*Wave Propagation and Scattering in Random Media*. New York: Academic Press, 572 p. Google Scholar

*Proc. SPIE 5885, Remote Sensing of the Coastal Oceanic Environment, 58850D*( San Diego, California). Kopilevich, Y.I.; Feygels, V.I.; Tuell, G.H., and Surkov, A., 2005. Measurement of ocean water optical properties and seafloor reflectance with scanning hydrographic operational airborne lidar survey (SHOALS): I. Theoretical background.

*Proc. SPIE 5885, Remote Sensing of the Coastal Oceanic Environment, 58850D*( San Diego, California).

*Journal of Optical Technology*, 75 (5), 321– 326. Google Scholar

*Journal of Optical Technology*, 77 (10), 598– 601. Google Scholar

*Proceedings of the VI International Conference Current Problems in Optics of Natural Waters*( St. Petersburg, Russia), pp. 183– 186. Kopilevich, Y.I.; Kononenko, M.E., and Zadorozhnaya, E.I., 2011. Prediction of bathymetric lidar performance with the account for the VSF shape.

*Proceedings of the VI International Conference Current Problems in Optics of Natural Waters*( St. Petersburg, Russia), pp. 183– 186.

*Proceedings of SPIE*, 7695, 76950Z-1– 12. Mathur, A.; Ramnath, V.; Feygels, V.; Fuchs, E.; Park, J.Y., and Tuell, G.H., 2010. Predicted lidar ranging accuracy for CZMIL.

*Proceedings of SPIE*, 7695, 76950Z-1– 12.

*In*: Bass, M. (ed.),

*Handbook of Optics, Second Edition*. New York: McGraw-Hill, 43.3– 43.56. Google Scholar

*Ocean optics. V. 1. Physical optics of the ocean*. Moscow: Nauka, in Russian. Google Scholar

*Proceedings of SPIE*, 7695, 769511, https://doi.org/10.1117/12.851978. Ramnath, V.; Feygels, V.; Kopilevich, Y.; Park, J.Y., and Tuell, G., 2010. Predicted bathymetric lidar performance of Coastal Zone Mapping and Imaging Lidar (CZMIL).

*Proceedings of SPIE*, 7695, 769511, https://doi.org/10.1117/12.851978.

*Marine Technology Society Journal*, 39 (3), 27– 35. Google Scholar