Derivation of the General Lidar Equation

Consider a lidar receiver at z₀-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:

Figure 2.

Mono-static lidar geometry, ‘R’ for receiver and ‘E’ for emitter.

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:

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

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:

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 dn_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:

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:

Calculation of the Normalized Irradiance Function

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

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

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₀ root

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⁻¹.

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²/4)P², 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).

Figure 3.

Characteristics of the components in the exponential term.

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).

Figure 4.

f(k;h,r) functions at h=5m, (a) Θ=2 mrad, (b) Θ=10 mrad.

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.

Figure 5.

Normalized irradiance functions.

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.

Figure 6.

Profiles of g_Θ (h,r) (a) for a Gaussian radiance source with varying distance h; (b) normalized at the peak; and (c) for Step-function source.

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.

Figure 7.

Irradiance dispersion (non-nadir incident mono-static system).

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 = (2r_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⁻¹(−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.

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 8.

(a) Normalized laser pulse. (b) Effect of system response.

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₀ 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⁻² 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.

Figure 9.

(a) CZMIL waveforms sampled from flat runway (circled area in the photo); (b) comparison with two simulated waveforms.

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.

Figure 10.

CZMIL radiometric calibration curve.

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.

Figure 11.

Waveforms from shallow to deep water along the direction marked in the image. Two deep-water bottom peaks are magnified.

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.

Figure 12.

CZMIL waveform and the simulated waveform for 2 m and 20 m.

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₁× n₂ + n₃ × n₄ ) / 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⁻¹(n_beam · n_bottom ) – θ_w ≈ 14.

Figure 13.

CZMIL waveform and the simulated waveform for 11 m.

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).

Figure 14.

CZMIL waveform and the simulated waveform for 30 m and 40 m.

Figure 15.

Analysis of F-function and system attenuation coefficients.

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.