1. Introduction

The understanding of the variability of optical properties in natural waters has important implications on mapping water quality parameters [1], productivity models [2], underwater communications [3], and underwater imaging [4]. Indeed, the spatial variability of inherent optical properties (IOPs) such as the beam attenuation coefficient (c) may be linked to turbidity plumes [5]. Likewise, photosynthesis and phytoplankton growth in aquatic environments is modulated by vertical light attenuation associated with c changes [6]. Lastly, c is a basic parameter for improving detection of submarine targets [7]. The remote sensing of c is only possible based on LiDAR (Light detection and ranging) systems [8]. Unlike passive optical sensors, LiDARs are able to measure characteristics of both the time-resolved back-scatter and forward-scatter photon contributions arriving at the detector. This is particularly feasible when the viewing angle and/or c magnitudes are relatively small [9,10]. In this case, the LiDAR waveform is mainly constituted by photons encompassing one collision over the return path. However, this is no longer true as water turbidity increases and/or the LiDAR field-of-view (FOV) becomes larger due to a greater contribution of photons going through multiple-scattering events.

To address these light propagation effects on IOPs inversion, different LiDAR equations have been proposed [9–12]. In that regard, there are two models commonly used in the literature: the quasi-single scattering [9,10] and the small-angle aproximation [11]. The first model was derived from the radiative transfer theory by assuming only one collision per photon. However, this premise is not strict as the attenuation term may also include a higher order of photon interactions. The second model is more suitable for relatively turbid waters and multiple scattering effects since it allows multiple forward scattering collisions at very small angles.

A common approach to extract c from waveforms obtained by oceanographic LiDARs is the calculation of one half of the slope of the log-transformed return as a function of range (hereafter α) [13]. Here, the factor of one half accounts for the two-way propagation of the return. In general, this mathematical procedure is straightforward and based on a linear regression fit of log-transformed volume backscattering values measured within the descending portion of the path-radiance peak. Despite its simplicity, this technique may result in ambiguous estimates if measurements are made across a wide range of water turbidities due to the variable and unknown contribution of multiple scattering. Thus, the 'slope inversion' may only be useful in clear waters and characterized by relatively small changes on IOPs.

An alternative way of retrieving c from LiDAR profiles is the characterization of the FOV loss function factor [7,14]. This parameterization describes the detector's efficiency for capturing photons with respect to c and range, and can be simulated or measured under controlled conditions [7]. Likewise, implicit solutions to the FOV loss function factor can be obtained by using multiple FOVs [15]. Therefore, the determination of the FOV loss function factor is not a straightforward task and can only be achieved by performing customized Monte Carlo-based radiative transfer simulations, time-consuming tank experiments, or having a complex LiDAR configuration with several FOVs.

Here, a new method to compute c from LiDAR waveforms produced by an underwater 3-D system is presented. The optical inversion is fast and based on an empirical model derived from a modified Weibull function that takes into account the time variation and signal/noise characteristics of each pulse return. Thus, the aim of this study is to quantify the response of Weibull parameters to c, and compare the resulting inversions with those derived from c-α relationships.

2. Data and methods

2.1 LiDAR system

The Fine Structure Underwater Imaging LiDAR (FSUIL) is a near-monostatic pulsed laser scanning system developed at Harbor Branch Oceanographic Institute [16]. The configuration of the orbital version includes a green (λ = 532 nm) and narrow collimated laser source with a divergence beam angle of 1 mrad, a 2-D scanner, and a receiver assembly with four 50 mm diameter telescopes (channel 1: wide-FOV and non-polarized, channel 2: wide-FOV and co-polarized, channel 3: wide-FOV and cross-polarized, channel 4: narrow-FOV and non-polarized). Each telescope has a bandpass filter centered at 532nm (3 nm at full width-half-maximum, FWHM, diameter = 50 mm), a F/2 plano-convex lens, a field stop iris, and a high speed photomultiplier tube (PMT) detector (Hamamatsu R9880U-210). The wide and narrow FOVs were controlled using the fixed iris at 75.7 and 15 mrad, respectively. The source-detector distance of channel 2 and 4 (0.157 m) is smaller with respect to that of channel 1 and 3 (0.266 m). Each pulse has a Gaussian spatial profile, a pulse duration at FWHM of 0.5 ns, and an energy of 20 µJ. The scanning area contains 460 pulses and covers an area of 1.7399 (width) x 1.5748 (height) m at 10.45 m from the detector. The scanning lines are created by swiping the laser beam from left to right for odd rows and viceversa at consecutive even rows. The angular inter-pulse separation is 0.34° and using a pulse repetition frequency (PRF) of 500Hz, the time to complete an individual scan is 0.92 seconds. The time to save the pulse matrix data to file is close to 2 seconds. Therefore the time between captures is approximately 3 seconds.

2.2 Experimental setup

FSUIL measurements were performed inside a water tank filled with a volume of freshwater of 151.42 cubic meters and specifically designed by Harbor Branch Oceanographic Institute for underwater optical experiments (Fig. 1). The tank’s dimensions are 6.5 x 12.5 x 2 m, and has two adjustable platforms on top to allow for precise deployment of instruments within the water volume. Extremely dark conditions are assured by having a completely light tight building, turning off artificial lights, and covering the surface of the water with black plastic spheres. The sphere cover and textured black walls and floor within the tank also work well in extinguishing internally scattered light, thereby creating a pseudo infinite volume for realistic ocean testing conditions. Before the experiments, 18 tank jets are directed downwards to resuspend any particles on the bottom. Likewise, potential residues resting at the bottom of the tank are eliminated by placing a series of filters in line (minimum pore size = 0.2 µm) and recirculating the water over a 24 hr period with a coagulation agent. The LiDAR measurements were done by placing the system in the center of the tank, near the edge and pointing along the major axis at a distant black target (2 m x 6 m, height x width) made of tightly stretched lint-free fabric, which was situated at 10.45 m with respect to the location of the FSUIL detector. A signal trigger is introduced to differentiate the beginning and the end of every waveform. The trigger is created by placing a beam sampler (Thorlabs) that splits the initial signal into two components. The small component (i.e., 1%) is going to the detector before even the transmitted pulse (i.e., the major component with 99% of the signal energy) leaves the housing of the LiDAR source.

 

Fig. 1 Test tank and LiDAR experimental settings.

Download Full Size | PDF

The influence of water turbidity on LiDAR waveforms was assessed based on 9 experiments where IOPs at λ = 532 nm were characterized by a mean value of c of 0.045, 0.236, 0.442, 0.660, 0.779, 0.980, 1.160, 1.330 and 1.520 m^-1, respectively. Values of a and c were measured using a WETLabs AC-9 instrument. Different water turbidities were obtained by adding a cumulative mass (i.e., 64 g per turbidity case study) of Ultra-fine Arizona Test Dust or ATD (Powder Technology, Inc.) to the volume of the tank. The dust was dropped from a bridge located near the target position. Notice that ATD in this study represents a model for suspended particulates. LiDAR measurements were initiated 15-25 minutes depending on final ATD concentration (i.e., longer at higher c values) and after adding the powder in order to reach a steady state ('mixed conditions') in terms of c homogeneity. The gain of each detector was 450 for channel 1 and 2, and 550 for channel 3 and 4.

2.3 Waveform simulations

The shape of the FSUIL waveforms was modeled empirically by implementing a modified expression of the Weibull probability distribution function (MW):

2.4 Statistical analysis

The Parameters P1, P2, P3, P4 and the LiDAR attenuation coefficient (α) were computed for each capture (i.e., single FSUIL scan) by choosing the pulse coinciding with the maximum intensity at the target location (i.e., pulse 288, 771, 1185 and 1644 for channel 1, 2, 3 and 4, respectively). This calculation was performed after smoothing of waveforms with a moving average of 10 samples. The magnitude of α was computed as half of the slope of napierian log-transformed and exponentially decreasing values of LiDAR backscattering as a function of range. The slope and its uncertainty were derived from type II linear regression models (i.e., independent and response variables are randomly selected). The experimental error of α and c estimates was quantified based on two metrics: the root mean square error (RMSE), and the median of unsigned relative differences (MURD):

For each α calculation, two points were dynamically selected from the trailing section of the path-radiance peak and based on the following criteria. First, the initial point or imax (i.e., value at an earlier time that corresponds to maximum) was found after eliminating the peaks corresponding to the trigger (TR) and target (TA), as shown in Fig. 2. Also, the code verifies for spikes using empirically-determined peak threshold, to avoid false detection of imax values. Spikes are detected if adjacent values to imax have a backscattering magnitude that is half of that observed at imax. Second, the final point or imin is computed as the first value above 1.5 or 3-fold the background signal (mean or maximum) of the last 15 time bins. Third, if the background signal is negative, then an offset is included as part of the imin detection. Lastly, if negative α values still occur then imin is computed as imax + 50 time bins, which was determined by visually inspecting the waveforms.

 

Fig. 2 Measured versus modeled FSUIL waveforms. a) channel 1, b) channel 2, c) channel 3, and d) channel 4. c is 0.045 m⁻¹ (red line), 0.236 m⁻¹ (black line) or 1.520 m⁻¹ (blue line); FSUIL measurements (solid lines), MW simulations (broken line), two standard errors (grey bars); TR and TA correspond to trigger and target positions, respectively.

Download Full Size | PDF

Similar to α, MW parameters were obtained after eliminating the trigger and target peaks. The uncertainty of Weibull function parameters was quantified based on two standard errors (i.e., standard deviation N^-0.5, where N is the number of observations) and the degree of function adjustment to the experimental data was evaluated based on ni or the number of iterations needed to reach the numerical convergence. The presence of outliers on calculating α, P1, P2, P3 and P4 was determined based on a Hampel filter [17].

3. Results

Unlike airborne LiDAR inversion algorithms typically applied in oceanographic applications, IOPs in this study are retrieved based on information derived from the whole waveform. In Fig. 2, examples of Weibull simulations are shown for low, intermediate and high concentrations of ATD across the four FSUIL channels. In this case, each LiDAR waveform corresponds to the arithmetic average of captures per experiment (i.e., 35, 25, and 44 scans for c = 0.045, 0.236 and 1.520 m⁻¹, respectively). Also for each curve, the position of TR and TA is indicated at 19 and 109 ns, respectively. Notice that backscattered intensity in channel 1 and 2 [Fig. 2(a) and 2(b)] was an order of magnitude greater than that measured in channel 3 and 4 [Fig. 2(c) and 2(d)]. These differences are explained by the smaller common scattering volume of channel 4 with respect to channel 1, 2 and 3. Also, a smaller contribution of path-radiance is expected in channel 1 with respect to channel 2 and 4 due to the longer source-detector distance of the former receiver.

In general for all channels, the leading edge of the path-radiance peak (i.e., first major bump before target return) was less variable and had a higher rate of change than the trailing edge. However, these differences were reduced at higher turbidities resulting in more symmetric Gaussian shapes. Also to notice was the presence of an additional path-radiance peak when the water turbidity was the lowest (i.e., c = 0.045 m⁻¹, red symbols in Fig. 2(a)). This feature showed earlier than the maximum path-radiance peak, was less remarkable in channel 3 and 4, and tended to disappear as c increases.

For all channels, with particular prominence in the lower attenuation cases, the peak associated with target-reflected photons is observed later than the path-radiance maximum and is still visible for the experiments made at c = 0.236 m⁻¹. Another clear feature in all waveforms except those measured at the highest turbidity, which were already significantly attenuated, was the drastic decay of the LiDAR backscattered power after hitting the black backwall target. Indeed, this signal variation was almost absent at c = 1.52 m⁻¹ due to the insignificant target return associated with a very strong path-length attenuation.

Regardless of the radiometric channel, the ratio between target and path-radiance photons (R_t-p) was only above 1 for the clearest waters. When visibility was the greatest, the maximum and minimum R_t-p values were obtained by using channel 4 (14.8) and 1 (1.02), respectively [Fig. 2(a) and 2(d)]. By keeping the FOV fixed, better discrimination of the target was achieved when a cross-polarized filter (i.e., channel 3) was used (R_t-p~7) [Fig. 2(c)].

A comparison of α calculations as a function of water turbidity and detector type is depicted in Fig. 3. As expected, the best linear adjustment between α and c values was found with measurements obtained by channel 4 [Fig. 3(d)]. The RMSE and MURD of c estimates for channel 1, 2, 3 and 4 was 0.147, 0.201, 0.253, and 0.127, and 13.7, 15.4, 15.3 and 9.4%, respectively. For the wide FOV, the correspondence between α and c values was greater when the unpolarized waveforms were analyzed [Fig. 3(a), 3(b) and 3(c)]. Indeed, the difference between coefficients of determination (r²) as percentage was up to 16% (Table 1).

 

Fig. 3 Variation of α as a function of c. a) channel 1, b) channel 2, c) channel 3, and d) channel 4. Linear regression model (broken lines), for each datapoint, vertical bars correspond to two standard errors.

Download Full Size | PDF

Table 1. Statistics of α-c relationships. Linear regression model: c’ = m₀ + m₁ α, c’ is modeled c, for each FSUIL channel, the number of observations was 9, two standard errors are indicated between parentheses

View Table | View all tables in this article

This was attributed to non-linearities of channel 2 at relatively high c values, and large variability around the regression line in channel 3. The range of α values for channel 1, 2, 3 and 4 was 0.081-0.186, 0.106-0.341, 0.091-0.202, and 0.056-0.320 m⁻¹, respectively. In general, the variability of α values among captures as inferred from the coefficient of variation (i.e., standard deviation/arithmetic average) as percentage was up to 14.6%, 28.8%, 28.8%, and 25.2% for channel 1, 2, 3 and 4, respectively, and tended to decrease at higher turbidities (up to 5.1% at c = 1.520 m⁻¹). In all cases, the error of α estimates as derived from FSUIL waveforms was always below 10%.

For all channels, the regression intercept of α as function of c was greater than zero (t- Student = 71.1, P < 0.01) and varied between 0.065 and 0.132. However, these differences were only apparent as they were comparable to the experimental error of c measurements (i.e., ~15%). Lastly, the regression slope of α as function of c and for channel 4 was 2.8-fold larger with respect to that computed for channel 1, 2 or 3 (Fig. 3).

The variation of FSUIL-derived Weibull parameters with respect to water turbidity is shown in Fig. 4. For the whole data set, the magnitude of P1 varied between 2.5 and 4.7, and had the largest uncertainty at relatively low water turbidities (i.e., c = 0.045-0.2 m⁻¹) when P1 decreased with c [Fig. 4(a)]. However, at higher water turbidities (i.e., c > 0.4 m⁻¹), P1 estimates had a linear and positive covariation with c. The arithmetic average of P1 per experiment was not substantially different among channels for c values between 0.2 and 1 m⁻¹. Nevertheless, a remarkable deviation characterized by smaller P1 at c > 1 m⁻¹ was found for channel 3 (Fig. 4(a), green symbols).

 

Fig. 4 Relationships between Weibull parameters and c. a) P1, b) P2, c) P3 and d) P4. Uncertainty bars correspond to two standard errors.

Download Full Size | PDF

The parameter P2 varied between 55.8 and 149.8, and presented a large variability among captures in relatively clear waters [Fig. 4(b)]. In general, the arithmetic average of P2 per experiment was larger for channel 1 and 3 with respect to channel 2 and 4. Despite these differences, all P2 estimates were related to c following a non-linear and inverse behavior. To model these functionalities, a third-order logarithmic function was proposed:

where c' is the c estimate and y_o, A, B and C are regression parameters summarized in Table 2 and computed based on the Marquardt-Levenberg algorithm.

Table 2. Statistics of P2-c relationships. Non-linear regression model: c' = y_o + A log_e P2 + B (log_e P2)² + C (log_e P2)³, c’ is modeled c, for each FSUIL channel, the number of observations was 9, two standard errors are indicated between parentheses

View Table | View all tables in this article

For all cases, this mathematical expression was able to explain more than 97% of the regression variability, and presented the largest accuracy (RMSE = 0.041, ni = 246) when waveforms were detected with channel 1 (Table 2).

The worst performance of Eq. (3) for predicting c was obtained with FSUIL measurements obtained with channel 2 (RMSE = 0.060, ni = 352).

As expected, all FSUIL waveforms were characterized by having larger path-radiance maximum values (i.e., P3) at higher water turbidities [Fig. 4(c)]. However, the rate of change of FSUIL intensity as a function of c, as inferred from the slope of the linear regression (m), was substantially greater for channel 2 (m = 2.33 ± 0.20, two standard errors) with respect to channel 1 (m = 1.30 ± 0.12), 3 (1.34 ± 0.40), and 4 (m = 1.07 ± 0.21). Another common pattern of all waveforms was the sigmoidal variation (i.e., a smaller rate of change at relatively low and high turbidities) of P3 with respect to c. This functionality was not well defined for channel 3 due to anomalous high backscattering at c = 0.236 m⁻¹ (see green symbols).

Noise of raw FSUIL waveforms (i.e., before moving averaging and P4 calculations) was not related to water turbidity and varied among channels with higher values in channel 1 (6.45 to 8.80) and 2 (0.71 to 3.92) with respect to channel 3 (−2.75 to −1.31) and 4 (−3.46 to −2.33). Due to the presence of negative intensity values, an offset equal to 5 was added to each waveform a priori of P4 inversions. This offset was only increased to 15 for channel 4 and the experiment was done with the clearest water.

4. Discussion

In this study, controlled experiments were performed to examine relationships between c and shape parameters derived from FSUIL waveforms. Strictly speaking, the validity of our results applied to waters with a c range of 0.045-1.52 m⁻¹. These water turbidity conditions are comparable to those used for measuring volume scattering functions (VSFs) in the Bahamas, Southern California coast, San Diego Harbor, Mediterranean Sea, Atlantic and Pacific Oceans, and Lake Baikal [18, 19]. Also, the suggested c-α and c-Weibull functionalities are probably correct in environments of the Tyrrenean Sea, and the English Channel since total scattering coefficient values reported by Kopelevich’s studies are consistent with those computed by Petzold’s VSF [20].

The use of an opaque target during the tank experiments is unnecessary for field measurements. However, this setup was used here in order to experimentally derive the FOV, and to find the pulse position having the maximum intensity in the far-field for each of the FSUIL channels. The presence of this target is not expected to affect the characterization of the path-radiance since its lower reflectivity (i.e., < 0.1) and depolarization properties [21]. Another aspect to highlight is the way turbidity case studies were created in the lab. Unlike FSUIL measurements in natural waters, tank experiments were done with suspended particulates having uniform properties, i.e. size distribution spectrum, refractive indices, and morphology. Conversely, IOPs distributions of field data are expected to have a shift toward large-sized particles as turbidity increases. Thus, the position of the maximum path-radiance observed here is likely to be displaced to a longer range (i.e., a later arrival time at the detector) if FSUIL is deployed in such a turbid natural environment. This effect can be attributed to the larger contribution of small-angle forward-scattered photons, the effect of which is an observable delay in the arrival of the backscattering path radiance. The presence of aggregates during the tank experiments cannot be ruled out and might explain the existence of spikes and anomalous values of Weibull inversions (e.g., P3 at c = 0.236 m⁻¹ in channel 3, Fig. 4(c)).

The analysis of α at different water turbidities clearly illustrated the suitability of a narrow-FOV for estimating c values, and the difficulty in applying a linear model to the co-polarized component of FSUIL measurements. Not surprising, the regression slope of c-α relationships had a maximum value in channel 4 as single-scattering effects dominate α. Since α is a variable constituted by c and K_d, and K_d or the vertically diffuse attenuation of downwelling irradiance is smaller than c, α functions depending on c are anticipated to have smaller regression slopes for wide-FOV channels (i.e., 1 to 3) as α tends to K_d due to a greater influence of multiple scattering [9].

As pointed by several studies [8, 9, 22], multiple-scattering effects are reduced when narrower FOVs such as the channel 4 are used. This effect explained the higher R_t-p values of channel 4 with respect to channel 1 in Fig. 2. This phenomenon is attributed to multiple photon collisions that enhance the path-radiance contribution and cause a reduction of the target component due to a stronger path-length attenuation of target-reflected photons. These changes are more remarkable in channel 1 with respect to channel 4.

The state of polarization of backscattered energy arriving at the FSUIL detector was also sensitive to changes on multiple-scattering due to variations on c. Indeed, the perpendicular component of polarization (i.e., channel 3) allowed larger R_t-p values with respect to the parallel component (i.e., channel 2). These differences were attributed to the larger contribution of multiple-scattering when LiDAR returns were co-polarized. The greater vertical penetration (i.e., smaller α) of LiDAR waveforms having a cross- with respect to a co-polarized component has been already reported for coastal waters [23]. In this case, LiDAR waveforms were obtained using an airborne system with a relatively narrow FOV (i.e., 17 mrad), and over surface waters encompassing c values between 1 and 5 m⁻¹ (Churnside, pers. Comm.).

Among all Weibull parameters under investigation, the parameter P2 was the most suitable for estimating c as a function of water turbidity due mainly to the simplicity of the inversion procedure, and weakly dependency of P2-based estimates with respect to the specific characteristics (e.g., transmitted power, dark current) of the LiDAR system.

The log-transformed polynomial in Eq. (3) had an outstanding predictive ability for deriving c and with respect to the traditional method based on the LiDAR backscattering slope. In fact, the best performance of P2-c parameterizations had an unexplained regression variability of less than 1% and much smaller with respect to that associated with c-α relationships (i.e., ~7%). Also, in terms of RMSE, P2-derived c estimates were up to 3-fold more accurate than those derived from α-based optical inversions. Although the superiority of using P2 for estimating c was clearly demonstrated, the exact mechanisms explaining c-P2 covariations were more elusive. It is suggested that as c increases, the contribution of path-radiance inside and outside the FOV increases due to multiple-scattering. Thus, as water became more turbid, the peaks corresponding to the first and the second path-radiance components merge resulting in a narrower path-radiance distribution (i.e., smaller P2 values).

The variation of P1 with respect to c had a dual trend showing an inverse and direct correspondence at relatively low (0.045-0.44 m⁻¹) and high (> 0.44 m⁻¹) c values, respectively. This complex pattern was attributed to shape modifications in the leading and trailing sections of the path-radiance maximum due to back-scattered and forward-scattered energy contributions that are originated inside (i.e., common volume backscattering) and outside of the FOV. In general, the proportion of photons having a multiple-scattering event is larger when the scattering process is initiated outside the FOV. Thus, it is suggested that P1 variations within the c range 0.045-0.236 m⁻¹ were mainly modulated by changes in path-radiance maximum (i.e., second path-radiance peak mainly composed of photons coming from the common volume backscattering). However at higher c values, the first path-radiance peak (i.e., ‘early-bird’ photons originated outside of FOV and close to the laser source) becomes more prominent and causes an augmentation of P1 (i.e., a symmetry improvement) as the first and maximum path-radiance peaks start merging. This explanation is supported by radiative transfer modeling and experimental results based on underwater laser line scanner measurements obtained with a LiDAR system having a comparable near-monostatic geometry to FSUIL [24].

5. Summary

In this investigation a novel and non-intrusive technique is proposed for estimating c based on LiDAR waveforms derived from FSUIL. The accuracy of P2-based models for estimating c was greater with respect to that typically obtained based on LiDAR backscattering slope inversions. Unlike α functions, the Weibull-based inversion utilizes the full waveform in order to calculate c. Also, another important benefit of using Weibull-derived retrievals was the fact that P2 can be used even if the waters under study are characterized by a variable multiple-scattering contribution. Future studies are underway for applying this model in oceanic and coastal waters. Likewise, the calculation of Weibull parameters for the entire set of pulses (i.e., 23 x 20) per capture will allow for the first time the fast reconstruction of 2-D and 3-D (if FSUIL is translated or rotated) distributions of c. These results are expected to have major applications in aquatic ecology (eg., patchiness) and defense (e.g., identification of submarine objects).

We thank HBOI engineering staff members Brian Ramos, Benjamin Metzger and Michael Young for electronic and mechanical design, fabrication and software development for the LiDAR instrument.