A Monte Carlo study of the seagrass-induced depth bias in bathymetric lidar

Chi-Kuei Wang; William Philpot; Minsu Kim; Hou-Meng Lei

doi:10.1364/OE.19.007230

1. Introduction

Airborne bathymetric lidar is a proven technology that can deliver bathymetric maps with high value to cost ratio. The system is mostly used in coastal zones where standard bathymetric mapping using conventional ship-based sonar is either too dangerous because of shallow depths and surf, or too expensive since multiple platforms (from small boats to large vessels) are needed to cover both the shallow and deep areas with any efficiency [1–3].

In determining depth from airborne lidar, it is generally assumed that the ocean bottom is clean and locally flat. In the natural environment, the ocean bottom varies significantly, both in depth and bottom characteristics, and the depth bias introduced by these factors is poorly understood. When seagrass is present, the greatest concern is whether the depth estimation will be biased [4,5]. The goal of this study is to investigate the effect of the presence of seagrass on the lidar waveform – the time history of the lidar return – from which water depth is estimated.

Lidar systems are widely employed both in ocean remote sensing and atmospheric studies [6,7]. Several analytical models for lidar systems are available, all of which rely on some degrees of simplification of the multiple scattering [7–11]. However, the complexity of representing the optical interactions on a seagrass bottom prohibits the implementation of any of these. A more practical approach is to use Monte Carlo simulation for the lidar signal [12–16].

2. Airborne bathymetric lidar system

For this study, the specifications of SHOALS1000 [17,18] are used to guide the modeling. The operational altitude range of SHOALS1000 is 200 - 400 m. The laser source transmits a 7 ns pulse at the repetition rate of 1 kHz. The energy scattered back from the water and the bottom can be stored at 1 ns resolution for data post processing. Table 1 lists the optical specifications of SHOALS1000 lidar system. In this study, the altitude of the lidar system is set at 400 m, at which SHOALS1000 usually operates. The transmitter and the receiver both are modeled as circular disks and located at the same position. The laser pulse enters the water at 20° nadir angle.

Table 1. 1 Optical Specifications of SHOALS1000 [18]

View Table | View all tables in this article

3. The Monte Carlo simulation

The essence of a Monte Carlo simulation is to follow a photon from the light source through all possible events (transmission, absorption, scattering, and reflection) that it may experience until reaching the receiver. This procedure is then repeated for each photon used for a simulation run. Each interaction event is determined by a properly associated probability function and the assigned physical properties of the materials (absorptivity, scattering phase function, etc.).

When a photon starts from the light source, it is assigned an initial weight of $w_{p}$ representing its initial energy. This quantity only represents the amount of energy carried by a virtual photon; it is not the energy of a photon in quantum physics. The distance traveled by the photon is recorded in order to allocate the photon energy in the appropriate time bin for lidar waveform construction. The lidar system digitizes the recorded signals into discrete intervals. These intervals are called time bins.

The Monte Carlo simulation for this study consists of several components to represent the natural environment and a specific lidar system: the photon, the optical hardware, the water body, the flat bottom, and the seagrass bed. The photon is traced through transmission and reflection at the air/water interface, scattering and absorption in the water, reflection and absorption at the ocean floor, and reflection, absorption, and transmission by the leaves of the seagrass. Atmospheric attenuation of the lidar signal is not considered because the operational altitude of the bathymetric lidar system is 400 m and the optical depth between the sensor and the ocean is assumed to be negligible.

The simulation begins with a photon emitted from the light source – in this case a laser – constrained by the transmitter FOV and the diameter of the light source. The details of the Monte Carlo technique used in this study are described below. A more detailed description of the technique can be found elsewhere [19,20].

In the description of the Monte Carlo processes, $w_{p}$ represents the photon energy before an interaction event with the air/water interface, water, seagrass leaves, or bottom while $w_{p}^{'}$ represents the photon energy after such an event.

3.1. Air/water interface processes

When the photon enters the air/water interface from above the water, the photon weight is reduced according to Fresnel reflectance, $r_{F}$

w_{p}^{'} = w_{p} (1 - r_{F}) .

And, its direction is changed according to Snell’s law of refraction.

When the photon contacts the air/water interface from below the water, three cases are tested. First, if the photon experiences total internal reflection, the photon weight does not change. Second, if the photon does not experience total internal reflection and would not be detectable by the sensor after exiting the water, only the reflected path is followed and the photon weight is reduced according to Fresnel reflectance. Last, if the photon will be detectable by the sensor after it exits the water, the photon weight is reduced by the Fresnel reflectance and is collected by the sensor and the photon is terminated. For the first two cases, the photon will change direction according to the law of reflection.

3.2. In water processes

The attenuation coefficient, c (m⁻¹), the scattering coefficient, b (m⁻¹), and the scattering phase function, $\tilde{β} (θ)$ (sr⁻¹), are the water IOPs (Inherent Optical Properties; the quantities that are independent of the illumination distribution), and θ is the scattering angle. After entering the water, the path length that the photon will travel is

l = - \ln (Λ),

where l is the optical path length (dimensionless) and Λ is a random number generated from a uniform distribution between 0 and 1. The optical path length is related to the physical path length, r (m⁻¹), such that l = cr. In this study, the water body is assumed to be vertically and horizontally homogenous. While this may not be the case for natural coastal ocean waters, it is a useful first step for this study since the primary purpose is to consider bottom effects, not water properties.

When the photon advances to the new location according to the randomly generated path length, if it is still in the water, i.e., no interactions with air/water interface or bottom materials (seagrass or ocean bottom), the photon energy is reduced by

w_{p}' = w_{p} \cdot ω_{0},

where

ω_{0}

is the single scattering albedo of water, defined as

ω_{0} = \frac{b}{c} .

The physical explanation of Eq. (3) is that the photon loses energy due to absorption by the water. At that point, the photon will be scattered with a reduced weight. The scattered direction of the photon is determined by its phase function $\tilde{β} (θ)$ , which describes the probability that a photon is scattered at angle θ. When integrated over all possible scattering directions the integral sums to unity:

2 π \int_{0}^{π} \tilde{β} (θ) \sin θ d θ = 1.

The water is modeled as consisting of two scattering components: pure sea water and particles. The phase function of pure sea water can be determined from Rayleigh theory:

{\tilde{β}}_{w} (θ) = \frac{3}{16 π} (1 + \cos^{2} θ) .

The particle phase function used here is that provided by [19], which is based on three phase function measurements of natural water from [21]. The phase function of the modeled sea water $\tilde{β} (θ)$ is then computed by combining a fraction of the phase function of pure sea water ${\tilde{β}}_{w} (θ)$ and that of particles ${\tilde{β}}_{p} (θ)$ .

\tilde{β} (θ) = \frac{{\tilde{β}}_{w} (θ) \cdot b_{w} + {\tilde{β}}_{p} (θ) \cdot b_{p}}{b_{w} + b_{p}},

where

b_{w}

and

b_{p}

are the scattering coefficient of pure sea water and particles, respectively.

Because Eq. (5) integrates to unity and has only one variable, the probability density function (pdf) of the scattering angle of a photon scattered from the water-particle mixture can be defined as

f_{θ} (θ) = 2 π \tilde{β} (θ) \sin θ, 0 \leq θ \leq π,

and the cumulative density function (cdf) can be defined as

F_{θ} (θ) = \int_{0}^{θ} 2 π \tilde{β} (θ^{'}) \sin θ^{'} d θ^{'}, 0 \leq θ \leq π .

The scattering angle of the scattered photon is determined using a new Λ and the relationship

θ = F_{θ}^{- 1} (Λ) .

In the simulation, the azimuth angle of the scattered photon is randomly chosen between 0 and 2π. So, another Λ is generated, and the azimuth angle of the scattered photon is 2πΛ.

After scattering in the water column or reflection from the bottom (see Section 3.3) photons must be traced back to the receiver. If each photon were followed until it reached the receiver, an unrealistic amount of computation time would be required. Because the receiver has a small aperture size and FOV, the probability of a single photon arriving at the receiver is almost zero. To improve the efficiency of the simulation, each photon in the water is assigned a fractional energy proportional to the probability that it will be scattered or reflected in the direction of the receiver [12,13,15,22]. For scattered photons we may write,

w_{w} = w_{P}^{'} \tilde{β} (θ_{T, W}) \times Δ Ω \times (1 - r_{F}^{'}) \times {(T_{G})}^{N} \times \exp (- c d_{T}),

where

θ_{T, W}

is the scattering angle from the photon direction

ξ^{'}

to the receiver, ΔΩ is the solid angle subtended by the receiver from the photon, d_T is the geometric distance in the sea water between the photon and the receiver,

r_{F}^{'}

is the Fresnel reflectance when the light is incident from water to air,

T_{G}

is the seagrass transmittance, and N is the number of seagrass leaves between the photon and the receiver. Figure 1(a) shows the geometric relationship of these quantities.

Fig. 1 The geometric relationship between the photon direction and the lidar receiver. The receiver collects a fraction of (a) scattered energy by sea water, (b) reflected energy by seagrass leaf, and (c) transmitted energy by seagrass leaf. [See the text for the explanation of the notations.]

Download Full Size | PDF

3.3. Ocean bottom processes

The ocean bottom is assumed to be a flat, Lambertian surface, which reflects the photon according to a cosine function. In the first attempt to investigate the influence of the presence of the seagrass, the bottom materials are modeled with this simplest assumption. A tilted bottom with a non-Lambertian surface can be added easily later on. Adding a bottom topology is possible, but will further complicate the model, would require more computational time, and is not necessary for the present purpose.

After reflection from the ocean bottom, the photon weight is reduced by,

w_{p}^{'} = w_{p} \cdot R,

where R is the bottom reflectance. The reflection angle

θ_{L}

of the reflected photon is defined as the angle between the normal vector of the bottom surface and the exit vector of the photon, i.e., 0° corresponds to exiting vertically and 90° corresponds to exiting horizontally. The cdf of

θ_{L}

is given by:

F_{L} (θ_{L}) = 2 \int_{0}^{θ_{L}} \cos θ^{'} \sin θ^{'} d θ^{'}, 0 \leq θ \leq \frac{π}{2} .

To determine a randomly chosen $θ_{L}$ , a new Λ is generated and $θ_{L}$ is determined as:

θ_{L} = F_{L}^{- 1} (Λ),

Another independent Λ is generated and the azimuth angle of the reflected photon is 2πΛ. As with the proportional energy collection of photons in the water, the weight contributed to the detected signal by a photon reflected from the bottom is

w_{B} = w_{P}^{'} \frac{\cos θ_{T, B}}{π} \times Δ Ω \times (1 - r_{F}^{'}) \times {(T_{G})}^{N} \times \exp (- c d_{T}),

where

θ_{T, B}

is the angle between the normal vector of the bottom surface and the direction from the photon to the receiver. The rest of the variables in Eq. (15) have the same definition as those in Eq. (11).

3.4. The seagrass bed processes

The seagrass canopy structure model is adopted from [23], which is deduced from field measurements of 32 turtlegrass (Thalassia testudinum Banks ex König) populations near Lee Stocking Island, Bahamas, and three eelgrass (Zostera marina L.) populations from California and Washington, USA. Both species of seagrass are modeled as elongated rectangular strips with a thickness of 250 μm and width of 5 cm [23]. In the Monte Carlo simulation, the seagrass is treated as if the blades had zero thickness in order to save the computation time while an absorption is assigned to the passage of the photon through the leaf. The temporal bias introduced by this simplification is on the order of one-millionth of a nano-second, which is much smaller than the sampling time of the normal lidar system, and is negligible. The nomenclature for describing seagrass structure is the same as in [23]. A sigmoid function is used to describe the vertical biomass distribution of the seagrass bed [23]:

B (z) = \frac{ψ}{1 + {[\frac{z}{I}]}^{s}},

where B(z) is the vertical biomass distribution (%), ψ is the basal biomass (%), z is the height above the seabed (m), I is the intermediate point within the canopy (m), and s is the shape factor for the canopy. Note that ψ is a quantity obtained by fitting the distribution of collected seagrass leaves using a mathematical model and extrapolating to z = 0 and I is a tuning variable for Eq. (16) without physical meaning. Basal biomass is related to the optically active biomass a few centimeters above seabed. However, it is strictly a fitting coefficient and cannot be measured directly. From the regression analysis of [23],

ψ = 2.51 h_{c}^{- 0.79}

and

I = 0.588 (1 - \exp (- 1.12 h_{c})),

where

h_{c}

is determined by the maximum canopy height (

h_{m}

) and the bending angle γof a seagrass bed by:

h_{c} = h_{m} \cos γ .

The bending angle of a seagrass leaf γ describes the tilt of a seagrass leaf away from the zenith, i.e., 0° corresponds to standing vertically and 90°corresponds to laying horizontally on the bottom. Since the seagrass leaf is assumed to consist of rectangular strips with homogeneous mass density, regardless of the different age stage of each leaf, Eq. (16) is essentially the pdf of the length of seagrass leaf in the bed.

It is worth mentioning that the biomass distribution is expressed as a percentage and is independent of the bending angles of the seagrass leaves. To determine the length of leaf according to the vertical seagrass biomass distribution and to avoid a complicated integration of Eq. (16), the acceptance-rejection method [24] is used, where two random numbers of uniform distribution between 0 and 1, Λ₁ and Λ₂, are generated independently. If

B (Λ_{1} \cdot h_{c}) \geq ψ \cdot Λ_{2},

the random number pair is accepted, and the leaf length is

Λ_{1} \cdot h_{c} / \cos γ

, otherwise, a new pair of random numbers is generated. The seagrass leaves are “planted” randomly in a 1 m x 1 m patch, which is then duplicated throughout the ocean bottom in order to conserve computer memory. For still water simulation, the azimuth angle of the seagrass leaf orientation is randomly chosen between 0 and 2π. To simulate the seagrass under the influence of the water current, all seagrass leaves are tilted to the same direction. The averaged value of 4.75 of the shape factor was obtained from the results of all seagrass beds in [23] and will be used for all the simulations for seagrass bed reconstruction as a representation of natural turtlegrass and eelgrass beds.

The reflection and transmission processes of the seagrass leaves are assumed to be Lambertian. Since no descriptions of the Bidirectional Reflectance Distribution Function (BRDF) or the Bidirectional Transmission Distributional Function (BTDF) of seagrass or monocot at leaf scale were found in the literature search, we justify the Lambertian assumption based on measurements and simulation of a dicot leaf. The major components of dicot leaf tissue consist of cuticle, epidermis, palisade cells, and spongy cells. The cuticle is the protective skin of the leaf and exists on both sides. Palisade cells are at the upper surface, while spongy cells are at the lower surface. Palisade cells serve as light guides to spongy cells. The two seagrasses under investigation are monocot leaves, which lack palisade cells.

In a ray-tracing model of a dicot leaf, with light incident at angle ranging from 0 to 60 degrees off nadir, the BTDF is close to that of a Lambertian transmitter [25]. In a reflectance measurement of the lower surface of a dicot leaf, with light incident at angles ranging from –30 to 30 degrees off nadir, the BRDF is close to that of a Lambertian reflector [26]. For both cases, spongy cells are at the exiting side of the leaf and act as an efficient light scatter. For seagrass leaves, full of spongy cells, the BRDF and BTDF can be assumed to be Lambertian.

In order to calculate the photon energy collected by the receiver, one must know whether the photon is reflected or transmitted to the receiver by the seagrass leaf. For the reflected case, the energy arriving at the receiver is

w_{G} = w_{P} \times R_{G} \frac{\cos θ_{T, G}}{π} \times Δ Ω \times (1 - r_{F}^{'}) \times {(T_{G})}^{N} \times \exp (- c d_{T}),

where

R_{G}

is the seagrass reflectance and

θ_{T, G}

is the angle between the normal vector of the seagrass leaf and the direction of the photon path to the receiver. For the transmitted case,

w_{G} = w_{P} \times T_{G} \frac{\cos θ_{T, G}}{π} \times Δ Ω \times (1 - r_{F}^{'}) \times {(T_{G})}^{N} \times \exp (- c d_{T}) .

Figures 1(b) and 1(c) show the geometric relationship of the quantities in Eq. (21) and Eq. (22). The variables in Eqs. (21) and (22) have the same definition as those in Eq. (11).

Before allowing the photon to proceed on to the next path segment, one must decide if the photon is reflected from or transmitted through the seagrass leaf. A new Λ is generated and compared with the ratio $ω_{G} = R_{G} / (R_{G} + T_{G})$ . The photon is reflected if Λ is less than $ω_{G}$ , otherwise it is transmitted though the leaf. The resultant photon weight for both cases is

w_{P}^{'} = w_{P} \times (1 - A_{G}),

where

A_{G}

is absorption of the seagrass leaf. The reflection/transmission angle of the photon is determined in a way similar to that of the reflection from the bottom.

3.5. Variance reduction and simulation verification

The result of the Monte Carlo simulation is the combination of all the possible events undergone by all of the photons that occur in the natural environment, each with a distinct probability. The final result, itself, has all the characteristics of a probability function. Hence, there is noise (variance) associated with the computed waveform. The resultant waveform will be more meaningful if the noise can be minimized. To reduce the noise, the simplest strategy is to increase the number of photons involved in the simulation. However, the noise decreases as a function of inverse square root of number of photons used [19]. For a time-invariant sensor, the noise will decrease to an acceptably low value with a reasonable increase of computation time. For the lidar system, however, the photon energy is distributed over the time span of the detector (1 ns intervals), which means that the probability of the photon arrives at the receiver within a particular time interval is quite low. A significant increase in the number of photons in the simulation – substantially increasing the required processing time – only results in a marginal decrease in noise.

A number of 10 million photons are used for each simulation through out this study. Two simulation results with different particle attenuation coefficients are shown in Fig. 2 (denoted as thin lines). It is a common feature in the simulated data that the basic lidar waveform seems to follow the lowest trend in the data and that most of the variance appears as positive spikes (higher return energy) as a result of photon counting process, which is essentially Poisson noise [27]. In order to remove the Poisson noise, we devise a median filter strategy, which conducts the simulation run nine times and the median operator is then applied to the nine resultant waveforms for each time bin. To further mitigate simulation noise and save computation time, each simulation is set up with the laser pulse length of 1 ns. The intermediate resultant median waveform is then convolved with the original pulse function (a 7 ns square function in this study) to produce the final waveform. Figure 2 shows the comparison between the normal simulation (photons fired within 7 ns pulse) and the simulation using the convolution of the median waveform (photons fired within 1 ns, repeated 9 times, apply the median operator to the 9 waveforms, the resultant median waveform is convolved with a 7 ns square pulse) for an ocean bottom at 20 m depth. The process is repeated for two different particle attenuation coefficients and the bottom peaks are collocated at the same position for both cases. The time starts from 0 when the 7 ns pulse starts firing.

Fig. 2 Comparisons of normal simulation (thin lines) and convolution of median waveform (thick lines) for beam attenuation coefficients of 0.25 and 0.45 m⁻¹. The bottom is 20 m deep.

Download Full Size | PDF

Furthermore, verification of the simulation results is achieved by comparing results of the simulation code with those for the PREMAR-2F model [22] using the same settings for the lidar system, water body, and ocean bottom (results not shown here). This confirms that the Monte Carlo code is suitable to investigating the bottom return of the lidar waveform.

4. Results and discussion

To investigate how the presence of seagrass affects the lidar waveform, the seagrass beds are laid out in two canopy height arrangements (0.5 and 1 m), three shoot densities (100, 500, and 1000 shoots/m²), and two leaf orientations with three different bending angles (5, 25, and 45 degrees). Two seagrass types, eelgrass and turtlegrass [23], and two bottom types, siliciclastic mud (referred as mud hereafter) [23] and carbonated sand (referred to as sand hereafter) [28] are the common substrates for eelgrass and turtlegrass, and are used for bottom materials. For seagrass in still water, the leaf azimuth angle is assigned randomly. For seagrass under the influence of a water current, the leaf azimuth angle is assigned with one single value. This results in a total of 288 simulations of seagrass beds. Figure 3 shows various examples of a 1 m × 1 m seagrass patch with 1 m canopy height and shape factor of 4.75 for two shoot densities, three bending angles and three water current conditions. The optical properties of seagrass and bottom materials are listed in Table 2 .

Fig. 3 Examples of seagrass beds (1 m x 1 m) with 1 m maximum canopy height and a shape factor of 4.75. Shoot density of 100 shoots/m² with bending angles of (A) 5° (B) 25° (C) 45°, when the seagrass is in still water. Shoot density of 1000 shoots/m² with bending angles of (D) 5° (E) 25° (F) 45°, when the seagrass is in still water. Shoot density of 1000 shoots/m² and azimuth angles of (G) 0° and (H) 90°, when the seagrass is under the influence of a water current. The arrows denote the moving direction of the simulated airborne lidar.

Download Full Size | PDF

Table 2. Optical Properties of Seagrass and Bottom Materials

View Table | View all tables in this article

As a demonstration of the effect of seagrass presence on the depth bias, one set of water IOPs is used for all of 288 simulation runs, where c = 0.25 m⁻¹, ω₀ = 0.6, and the ocean bottom is set at the depth of 9 m.

The canopy density presented to the lidar varies with the actual shoot density, the bending angle of the leaves, and the orientation of the leaves relative to the lidar. The standard estimate of vegetation density is the leaf area index (LAI), the ratio of the projected leaf area on the bottom and the bottom area. LAI is useful when the sensor is viewing the seagrass from the zenith angle. The lidar receiver is viewing the bottom at a slant angle. Instead, an effective LAI (eLAI) is used. The eLAI is defined as the ratio of the projected leaf area on the plane perpendicular to the lidar incident vector and the bottom area projected to that plane.

The same seagrass structure, hence, the same LAI value, will have different eLAI values for different leaf azimuth angles. It can be defined as:

eLAI = \frac{A_{G} \times \cos θ_{i, g}}{A_{B} \times \cos θ_{i}},

where

A_{G}

is the total leaf area within unit patch,

A_{B}

is the area of unit patch,

θ_{i, g}

is the incident angle on the seagrass leaf (relative to the normal of the seagrass leaf), and

θ_{i}

is the refracted angle of the laser beam for a flat ocean bottom with zero slope. The eLAI is computed for the simulated seagrass bed for each simulation run.

It should be noted that the length of each seagrass leaf (c.f. Fig. 3) is determined based on the mathematical model being used (Eq. (16)), which reflects the leaf structure of a seagrass canopy. Due to the low probability of producing a leaf that reaches the maximum height specified for the canopy, the longest leaves generated by the simulation are sometimes shorter than the maximum height specified by the user.

The depth estimation algorithm of SHOALS1000 determines the time bin of the half peak value of the bottom return by subtracting the estimated water return from the lidar waveform [3]. To avoid developing our own signal-processing algorithm, we take advantage of having full control of the simulation. All the received photons that have been reflected from the ocean bottom or have interacted with seagrass are grouped into the bottom return of the simulation, and the time bin of the half peak value of the simulated bottom return is used for depth estimation.

The depth bias $Δ D$ (ns) of each simulation run is quantified by the difference between the half peak time bin of the bottom return, $t_{G}$ (ns), and that of a flat bottom without seagrass $t_{B}$ (ns):

Δ D = t_{G} - t_{B} .

The summarized results, including the maximum, minimum and mean eLAI values, for different seagrass parameters are listed in Table 3 . For all of the simulation runs, the maximum and minimum eLAI are 1.56 and 0.01, respectively, and the eLAI of turtlegrass is similar to that of eelgrass. The eLAI generally increases with increasing shoot density, bending angle, and canopy height. The seagrass has the largest and smallest eLAI values when the leaf azimuth angles are 0 (the seagrass leaves are facing the laser) and 180 degrees, respectively. The eLAI of seagrass in still water (random leaf azimuth angle) is similar to that for oriented leaves with an azimuth angle of 90 degrees, which implies that the seagrass with randomly distributed leaf azimuth angles can be presented as a mean leaf azimuth angle of 90 degrees. Note that, the range (from minimum to maximum value) of eLAI overlap within each seagrass parameter. This is due to the fact that there are other seagrass structure variations within each seagrass structure setting. For example, the maximum eLAI of seagrass with bending angle of 5 degrees (shoot density of 1000 shoots/m², canopy height of 1 m, leaf azimuth angle of 0 degree) is larger than the minimum eLAI of that with 45 degrees (shoot density of 100 shoots/m², canopy height of 0.5 m, leaf azimuth angle of 180 degrees).

Table 3. Summary of eLAI and Depth Error for Different Seagrass Parameters

View Table | View all tables in this article

To illustrate the rusults, a subset of the 288 simulation runs is shown in Fig. 4 , which shows results for turtlegrass planted on sand and mud bottoms with leaf azimuth angle of 0 degree, canopy height of 1 m, with the bending angle varying from 5 to 45 degrees, and shoot density varying from 100 to 1000 shoots/m². In general, other seagrass parameters being the same, the depth bias increases with shoot density and bending angle (Table 3). However, one single seagrass parameter is insufficient to represent the depth bias. A better correlation relationship is found between depth bias and eLAI, where the depth bias increases negatively (meaning the bottom determined by lidar is shallower than the bottom without seagrass) with increasing eLAI. Two trends of depth bias for sand (open symbol) and mud (solid symbol) bottoms can be seen in Fig. 4. It is postulated that the reflectance contrast between sand bottom and turtlegrass is not as significant as that between mud bottom and turtlegrass (e.g., Table 2), so the depth bias induced by similar eLAI values (i.e., similar seagrass structural and optical properties) is less for a sand bottom. Thus, the summarized results of the depth bias for different seagrass parameters listed in Table 3 are grouped by bottom type.

Fig. 4 A subset of the results showing the depth error induced by turtlegrass with leaf azimuth angle 0° and canopy height of 1 m. The open and solid symbols denote sand and mud bottoms, respectively; the triangles, circles, and squares denote the bending angles of 5°, 25°, and 45°, respectively; red, blue, and green denote the shoot density of 100, 500, and 1000 shoots/m², respectively.

Download Full Size | PDF

When eLAI is smaller than 0.07, the maximum values of the depth bias for sand and mud bottoms and for various seagrass parameters is very close to zero, meaning that the depth bias is negligible (Table 3). Nonetheless, the depth bias is more pronounced for the mud bottom than for the sand bottom (the minimum values of depth error of mud bottom are smaller than those of sand bottom). Based on the minimum and mean values of the depth bias in Table 3, for both bottom types, the depth bias induced by turtlegrass and eelgrass is similar; the depth bias increases negatively with increasing shoot density, canopy height, and bending angle.

The depth biases are most and least pronounced when leaf azimuth angles are 0 and 180 degrees, respectively. Also, the depth bias for random leaf azimuth angles is similar to that of leaf azimuth angle of 90 degrees, which is not surprising because their eLAI values are quite similar.

Figure 5 shows the results of all 288 simulation runs. For ease of comprehension, Figs. 5(a) and 5(b) are color- and symbol- coded differently according to the optical parameters in Table 2.

Fig. 5 The depth errors induced by seagrass. (a) and (b) shows exactly the same data but are color- and symbol- coded differently according to parameters setup in Table 2. (a) The open and solid symbols denote sand and mud bottoms, respectively; the squares and circles denote eelgrass and turtlegrass, respectively; red denotes the seagrass in still water and green, blue, and cyan denote the seagrass affected by water current of 0°, 90°, and 180°, respectively. (b) The open and solid symbols denote the seagrass height of 0.5 m and 1 m, respectively; the triangles, circles, and rectangles denote the bending angle of 5°, 25°, and 45°, respectively; the red, blue, and green denote shoot density of 100, 500, and 1000 shoots/m², respectively.

Download Full Size | PDF

The depth bias induced by turtlegrass and eelgrass is similar (shown as pairs of circles and squares in Fig. 5(a)) because the optical properties of turtlegrass are similar to those of eelgrass. The two trends of depth biases for sand (open symbols) and mud (solid symbols), shown in Fig. 4, are also seen in Fig. 5(a). If the depth bias of −1 ns is arbitrarily used as a threshold, only one case for sand bottom (shoot density of 1000 shoots/m², canopy height of 1 m, bending angle of 45 degrees, and leaf azimuth angle of 0 degree; shown as green open circle in Fig. 5(a) and green solid rectangle in Fig. 5(b)) would be considered to be significantly biased. On the other hand, many of the cases for mud bottom have the depth bias smaller than −1 ns, and some are biased as much as −2.22 ns.

To further utilize the results shown in Fig. 5, we suggest data collection of in situ measurements, concurrent with the airborne bathymetric lidar survey, including water IOPs, bottom (substrate) reflectance, and seagrass structure parameters (i.e., seagrass type, shoot density, maximum canopy height, bending angle and leaf azimuth angle) [23] for depth correction purpose.

5. Conclusions

A Monte Carlo simulation including a plane parallel water body, a flat bottom, and seagrass, is demonstrated in this study. The seagrass canopy structure model is adopted from [23], which is deduced from field measurements of 32 turtlegrass (Thalassia testudinum Banks ex König) populations near Lee Stocking Island, Bahamas, and three eelgrasses (Zostera marina L.) populations from California and Washington, USA. Both seagrasses are modeled as elongated rectangular strips with a thickness of 250 μm and width of 5 cm [23]. The bottom and seagrass leaf are considered lambertian surface. To save computation time, a median filter strategy accompanying with waveform convolution is adopted.

For all 288 simulations runs considered in this study, it is found that the depth bias is correlated with seagrass structure and reflectance contrast between seagrass and the bottom. Furthermore there is no depth bias difference between turtlegrass and eelgrass. By characterizing seagrass structure using eLAI, a more pronounced correlation relationship is found between eLAI and depth bias (Figs. 4 and 5). Because the reflectance contrast is more pronounced for seagrass and mud bottom, more significant depth bias is found for mud bottom than sand bottom.

Acknowledgments

Chi-Kuei Wang is supported by the National Science Council of Taiwan (Grant Number NSC 99-2628-E-006-160). We thank the anonymous referee for the constructive comments.

References and links

1. G. C. Guenther, “Airborne lidar bathymetry,” in Digital Elevation Model Technologies and Applications: The DEM Users Manual, D. F. Maune, ed. (ASPRS, 2001)

2. G. C. Guenther, A. G. Cunningham, P. E. LaRocque, and D. J. Reid, “Meeting the accuracy challenge in airborne lidar bathymetry,” EARSeL eProceedings (2001), Vol. 1, pp. 1–27.

3. G. C. Guenther, R. W. L. Thomas, and P. E. LaRocque, “Design considerations for achieving high accuracy with the shoals bathymetric lidar system,” Proc. SPIE 2964, 54–71 (1996). [CrossRef]

4. M. Lee, “Benthic mapping of coastal waters using data fusion of hyperspectral imagery and airborne laser bathymetry,” PhD dissertation (University of Florida, 2003).

5. C.-K. Wang and W. D. Philpot, “Using airborne bathymetric lidar to detect bottom type variation in shallow waters,” Remote Sens. Environ. 106(1), 123–135 (2007). [CrossRef]

6. L. R. Bissonnette, G. Roy, L. Poutier, S. G. Cober, and G. A. Isaac, “Multiple-scattering lidar retrieval method: tests on Monte Carlo simulations and comparisons with in situ measurements,” Appl. Opt. 41(30), 6307–6324 (2002). [CrossRef] [PubMed]

7. R. M. Measures, Laser Remote Sensing (Krieger, 1992)

8. L. R. Bissonnette, “Multiple scattering of narrow light beams in aerosols,” Appl. Phys. B 60(4), 315–323 (1995). [CrossRef]

9. C. Flesia and P. Schwendimann, “Analytical multiple-scattering extension of the Mie theory: the LIDAR equation,” Appl. Phys. B 60(4), 331–334 (1995). [CrossRef]

10. A. V. Starkov, M. Noormohammadian, and U. G. Oppel, “A stochastic model and a variance-reduction Monte-Carlo method for the calculation of light transport,” Appl. Phys. B 60(4), 335–340 (1995). [CrossRef]

11. E. P. Zege, I. L. Katsev, and I. N. Polonsky, “Analytical solution to LIDAR return signals from clouds with regard to multiple scattering,” Appl. Phys. B 60(4), 345–353 (1995). [CrossRef]

12. P. Bruscaglioni, A. Ismaelli, and G. Zaccanti, “Monte-Carlo calculations of LIDAR returns: procedures and results,” Appl. Phys. B 60(4), 325–329 (1995). [CrossRef]

13. H. R. Gordon, “Interpretation of airborne oceanic lidar: effects of multiple scattering,” Appl. Opt. 21(16), 2996–3001 (1982). [CrossRef] [PubMed]

14. G. W. Kattawar and G. N. Plass, “Time of flight lidar measurements as an ocean probe,” Appl. Opt. 11(3), 662–666 (1972). [CrossRef] [PubMed]

15. L. R. Poole, “Radiative transfer model for airborne laser fluorosensors: inclusion of water Raman scattering,” Appl. Opt. 21(17), 3063–3065 (1982). [CrossRef] [PubMed]

16. D. M. Winker and L. R. Poole, “Monte-Carlo calculations of cloud returns for ground-based and space based LIDARS,” Appl. Phys. B 60(4), 341–344 (1995). [CrossRef]

17. J. Heslin, W. J. Lillycrop, and R. Pope, “CHARTS: an evolution in airborne lidar hydrography,” presented at U.S. Hydro Conference, Biloxi, Missippi, 24–27 March 2003.

18. G. Cunningham, Marine Survey Division, Optech Inc., 100 Wildcat Road, Toronto, Ontario M3J 2Z9, Canada (personal communication, 2004).

19. C. D. Mobley, Light and Water: Radiative Transfer in Natural Waters (Academic, 1994).

20. G. I. Marchuk, G. A. Mikhailov, M. A. Nazaraliev, R. A. Darbinjan, B. A. Kargin, and B. S. Elepov, The Monte Carlo Methods in Atmospheric Optics, (Springer-Verlag, 1980)

21. T. J. Petzold, “Volume scattering functions for selected ocean waters,” SIO Ref. 72–78 (Scripps Institution of Oceanography, 1972).

22. R. Barbini, F. Colao, E. Cupini, N. Ferrari, G. Ferro, and A. Palucci, “Marine code for modelling range resolved oceanographic lidar fluorosensor measurements,” EARSeL eProceedings (2001), Vol. 1, pp. 77–87.

23. R. C. Zimmerman, “A biooptical model of irradiance distribution and photosynthesis in seagrass canopies,” Limnol. Oceanogr. 48(1_part_2), 568–585 (2003). [CrossRef]

24. R. Y. Rubinstein, Simulation and the Monte Carlo Method (Wiley, 1981)

25. Y. M. Govaerts, S. Jacquemoud, M. M. Verstraete, and S. L. Ustin, “Three-dimensional radiation transfer modeling in a dicotyledon leaf,” Appl. Opt. 35(33), 6585–6598 (1996). [CrossRef] [PubMed]

26. T. W. Brakke, J. A. Smith, and J. M. Harnden, “Bidirectional scattering of light from tree leaves,” Remote Sens. Environ. 29(2), 175–183 (1989). [CrossRef]

27. R. Zanella, P. Boccacci, L. Zanni, and M. Bertero, “Efficient gradient projection methods for edge-preserving removal of Poisson noise,” Inverse Probl. 25(4), 045010 (2009), doi:. [CrossRef]

28. E. J. Hochberg, M. J. Atkinson, and S. Andrefouet, “Spectral reflectance of coral reef bottom-types worldwide and implications for coral reef remote sensing,” Remote Sens. Environ. 85(2), 159–173 (2003). [CrossRef]

Transmitter – 532 nm		Receiver – 532 nm
Beam divergence	5 mrad @ 400 m	FOV (Field of View)	21 mrad
Pulse width	7 ns	Aperture size	100 mm
Output power	5 mJoule

Seagrass type	Reflectance (R_G)		Transmittance (T_G)
Turtlegrass^a	0.0607		0.318
Eelgrass^a	0.0659		0.281
Seagrass parameter		Values
Bending angle [°]		5, 25, 45
Maximum canopy height [m]		0.5, 1
Shoot density [shoots/m²]		100, 500, 1000
Leaf azimuth angle [°]		0, 90, 180, random
Bottom type	Reflectance (R)
Sand^b	0.51
Mud^a	0.0364

		eLAI			Depth Bias [ns]
		eLAI			Sand Bottom			Mud Bottom
		max	min	mean	max	min	mean	max	min	mean
Seagrass type	Eelgrass	1.56	0.01	0.29	0.05	−0.97	−0.05	0.00	−2.22	−0.38
Seagrass type	Turtlegrass	1.51	0.01	0.29	0.07	−1.02	−0.06	0.00	−2.19	−0.41
Shoot density [shoots/m²]	100	0.17	0.01	0.06	0.04	−0.05	0.00	0.00	−0.40	−0.08
	500	0.76	0.03	0.27	0.05	−0.29	−0.03	−0.01	−1.46	−0.38
	1000	1.56	0.06	0.54	0.07	−1.02	−0.14	−0.03	−2.22	−0.72
Maximum canopy height [m]	0.5	0.77	0.01	0.22	0.04	−0.14	−0.01	0.00	−0.84	−0.21
Maximum canopy height [m]	1	1.56	0.02	0.35	0.07	−1.02	−0.10	0.00	−2.22	−0.58
Bending angle [°]	5	0.27	0.01	0.10	0.04	−0.08	−0.01	0.00	−0.32	−0.10
	25	0.90	0.02	0.24	0.07	−0.36	−0.01	0.00	−1.52	−0.31
	45	1.56	0.07	0.52	0.02	−1.02	−0.14	−0.05	−2.22	−0.78
Leaf azimuth angle [°]	random	1.22	0.01	0.29	0.02	−0.62	−0.05	0.00	−2.00	−0.39
	0	1.56	0.02	0.35	0.03	−1.02	−0.10	0.00	−2.22	−0.50
	90	1.23	0.01	0.28	0.03	−0.58	−0.05	0.00	−2.15	−0.40
	180	0.92	0.01	0.23	0.07	−0.16	−0.01	0.00	−1.60	−0.29

Transmitter – 532 nm		Receiver – 532 nm
Beam divergence	5 mrad @ 400 m	FOV (Field of View)	21 mrad
Pulse width	7 ns	Aperture size	100 mm
Output power	5 mJoule

Seagrass type	Reflectance (R_G)		Transmittance (T_G)
Turtlegrass^a	0.0607		0.318
Eelgrass^a	0.0659		0.281
Seagrass parameter		Values
Bending angle [°]		5, 25, 45
Maximum canopy height [m]		0.5, 1
Shoot density [shoots/m²]		100, 500, 1000
Leaf azimuth angle [°]		0, 90, 180, random
Bottom type	Reflectance (R)
Sand^b	0.51
Mud^a	0.0364

A Monte Carlo study of the seagrass-induced depth bias in bathymetric lidar

Abstract

1. Introduction

2. Airborne bathymetric lidar system

3. The Monte Carlo simulation

3.1. Air/water interface processes

3.2. In water processes

3.3. Ocean bottom processes

3.4. The seagrass bed processes

3.5. Variance reduction and simulation verification

4. Results and discussion

5. Conclusions

Acknowledgments

References and links

Cited By

Figures (5)

Tables (3)

Equations (25)

Optics Express