Interpretation of Spectral LiDAR Backscattering off the Florida Coast

: A multispectral backscattering LiDAR (Light detection and range) system (hereafter Oculus) was integrated into a wave glider and used to estimate the scattering order (i.e., single vs multiple collisions) of LIDAR backscattering, the water inherent optical properties (IOPs), the biogeo-chemical characteristics of particulate scatterers (i.e., relative size, composition) and their motion) on shelf waters of South East Florida. Oculus has a dual-wavelength conﬁguration (473 and 532 nm) and two detection geometries (off- and on-axis). Characteristics of scatterers were investigated based on two complementary LiDAR-derived proxies (the Structural Dissimilarity Index and the spectral slope of LiDAR backscattering). In March 2017, ﬁeld measurements showed a covariation between direct and diffuse backscattering contributions during morning hours and away from shore. LiDAR attenuation coefﬁcients explained up to 57% of IOPs variability. The analysis of LiDAR-derived proxies suggested higher turbidity and larger particulates near the coast


Introduction
The characterization of underwater scatterers based on light and range detection (LiDAR) measurements has been fundamental in studies related to the mapping of turbidity plumes [1] and thin scattering layers [2]. The main finding of these contributions was the differentiation of scattering layers in terms of vertical (e.g., nepheloid vs water column) and horizontal distributions (e.g., plankton patchiness). In that regard, the composition of scattering layers has been largely unknown for more than one decade due in part to the poor spectral resolution of LiDAR systems for water applications. To cope with this limitation, different techniques based on hybrid information (e.g., spectral reflectance and LiDAR backscattering) [3], relationships between optical properties derived from LiDAR waveforms [4], spatial statistics of LiDAR backscattering magnitude [5], signal thresholds (e.g., detection of fish schools) [6] and complementary use of hydrodynamic model simulations (e.g., Langmuir cells) [7] have been reported.
The accurate detection and identification of relatively large scatterers (i.e., size parameter = π D/λ >> 1, where D is the scatterer diameter and λ is the wavelength) [8] highly relies on how well the 'background' scattering of the optical medium is removed. This baseline signal is determined by the inherent optical properties (IOPs) of the waters under investigation and is critical for optimizing LiDAR-based imagers of underwater features [9]. Also, LiDAR-derived IOPs (e.g., beam attenuation coefficient, c) can be used in models for estimating signal scattering orders (i.e., single vs multiple collisions) [10] and apparent optical properties (e.g., K d or diffuse attenuation coefficient of downwelling irradiance, Table 1) needed on image denoising [4]. The use of autonomous robotic platforms such as gliders has major potential for mapping large-scale (i.e., 500 km) distributions of optical properties (e.g., total backscattering coefficient (b b )) and derived biogeo-optical variables (e.g., chlorophyll-a concentration, phytoplankton composition) in marine waters at high spatial resolution (i.e., cm) [11]. Despite their scientific value, these glider-based optical determinations may be influenced by the light field (e.g., upward measurements near the surface during daytime) and gliderassociated turbulence. Likewise, existing glider-based optical sensors do not provide range-resolved information and are unable to distinguish relatively large scatterers.
Here, a multi-FOV (Field-Of-View) and dual-wavelength LiDAR system (hereafter Oculus, Figure 1) is evaluated for characterizing IOPs and large scatterers in shelf waters of SE Florida. Oculus was developed for NOAA (National Oceanic and Atmospheric Administration) for studying the behavior of marine life and can be deployed in wave gliders (i.e., wave-propelled robotic platforms) [12]. The main advantages of using Oculus for detecting and discriminating scatterers are high spatial/temporal resolution (i.e., 100 waveforms per s, vertical resolution = 5.625 cm) [10] and the existence of two receivers (on-and off-axis) for estimating direct and diffuse scattering contributions to the total backscattering signal.
This study has three main goals: (1) to determine diel changes of direct and diffuse scattering components with respect to the distance to the shore, (2) to examine relationships between LiDAR optical properties (e.g., LiDAR or system attenuation coefficient, K sys ) and IOPs for shelf waters of South Florida having a wide range of turbidity (i.e., range of c at a wavelength of 532 nm 0.02-0.5 m −1 ), and (3) to evaluate two LiDAR-based and complementary proxies (the Structural Dissimilarity Index, SDI [13] and the spectral slope of LiDAR backscattering, m k [14] for discriminating different scatterers in terms of motion, relative size and composition. This contribution is organized into three main sections. In Section 1, LiDAR scattering contributions (direct and diffuse) are estimated based on Independent Component Analysis (ICA) [15] of waveforms arriving to off-and on-axis receivers. Also, the variability of these contributions in our study area was interpreted with respect to environmental factors. In Section 2, Oculus-derived K sys values were related to c and total absorption coefficients (a) in relatively clear and turbid waters (i.e.,~c H = 3 optical depths, where H is the water depth in m). Lastly in Section 3, scatterers were classified based on 2-D structure patterns (spatial and temporal) and spectral backscattering changes linked to particle size spectra and organic/inorganic content of suspended particulates.

The LiDAR System
The instrument has two non-scanning lasers with wavelengths centered in the blue (λ = 473 nm) and green (λ = 532 nm) spectral range ( Figure 1) [10]. The laser beams are parallel and were oriented at a nadir angle (i.e., perpendicular to the surface of the water and down-looking) ( Figure 2).
The receivers have an angular separation of 174.5 mrad and a similar FOV of 34.9 mrad. This geometry enhances the detection and discrimination of optical targets by allowing a real-time baseline correction by measuring direct and diffuse backscattered photons (i.e., path-radiance) in a concurrent way. The beam divergence is 17.5 mrad and the source-receiver separation is 0.0606 m for both telescopes. This source-receiver (S-R) geometry is an important design feature as high frequencies (i.e., >10 8 Hz) and associated dephasing of backscattered photons are sensitive to changes in S-R distance [16]. Oculus has an averaged laser power of 10 mW and 23.1 mW at 473 nm and 532 nm, respectively, a pulse repetition rate of 100 Hz and a pulse length of 1.27 and 1.12 ns (blue and green channels, respectively). The sampling frequency (i.e., digitization rate) during all experiments was 0.5 GHz.

The LiDAR System
The instrument has two non-scanning lasers with wavelengths centered in the blue (λ = 473 nm) and green (λ = 532 nm) spectral range ( Figure 1) [10]. The laser beams are parallel and were oriented at a nadir angle (i.e., perpendicular to the surface of the water and down-looking) ( Figure 2).
The receivers have an angular separation of 174.5 mrad and a similar FOV of 34.9 mrad. This geometry enhances the detection and discrimination of optical targets by allowing a real-time baseline correction by measuring direct and diffuse backscattered photons (i.e., path-radiance) in a concurrent way. The beam divergence is 17.5 mrad and the sourcereceiver separation is 0.0606 m for both telescopes. This source-receiver (S-R) geometry is an important design feature as high frequencies (i.e., >10 8 Hz) and associated dephasing of backscattered photons are sensitive to changes in S-R distance [16]. Oculus has an averaged laser power of 10 mW and 23.1 mW at 473 nm and 532 nm, respectively, a pulse repetition rate of 100 Hz and a pulse length of 1.27 and 1.12 ns (blue and green channels, respectively). The sampling frequency (i.e., digitization rate) during all experiments was 0.5 GHz.

Field Experiments
Optical measurements were made on 10 March 2017 during daytime conditions. Surveys were done over shelf waters off Fort Lauderlade, Florida (26.1224 • N, −80.1373 • W). LiDAR measurements were obtained from a rotating pole attached to the stern of the ship (Newton 40) ( Figure 2). LiDAR surveys were complemented with vertical profiles of c and a coefficients as derived from an absorption-attenuation meter (ac-9, accuracy ±0.001 m −1 , sampling rate = 3 Hz, SEA-BIRD Scientific, Bellevue, WA, USA) at 9 wavelengths (λ = 412, 443, 488, 510, 532, 555, 650, 676 and 715 nm). Protocols describing processing and signal corrections applied to raw ac-9 measurements are reported in previous studies [17].
The sampling design encompassed a time series over relatively deep waters (i.e., bottom depth range = 13.

Field Experiments
Optical measurements were made on 10 March 2017 during daytime conditions. Surveys were done over shelf waters off Fort Lauderlade, Florida (26.1224°N, −80.1373°W). LiDAR measurements were obtained from a rotating pole attached to the stern of the ship (Newton 40) ( Figure 2). LiDAR surveys were complemented with vertical profiles of c and a coefficients as derived from an absorption-attenuation meter (ac-9, accuracy ±0.001 m −1 , sampling rate = 3 Hz, SEA-BIRD Scientific, Bellevue, WA, USA) at 9 wavelengths (λ = 412, 443, 488, 510, 532, 555, 650, 676 and 715 nm). Protocols describing processing and signal corrections applied to raw ac-9 measurements are reported in previous studies [17].

Direct and Diffuse Scattering Components
The backscattering signals arriving to on-and off-axis receivers can be interpreted as a linear sum of two variable backscattering contributions (direct and diffuse) associated with single and multiple scattering collisions, respectively. In that regard, ICA deals with the partition of a mixed-signal that is composed of linear contributions terms. Unlike principal components, ICA is assuming latent variables having 'non-Gaussian' distributions and independent components that are not necessarily orthogonal. Another important as-

Direct and Diffuse Scattering Components
The backscattering signals arriving to on-and off-axis receivers can be interpreted as a linear sum of two variable backscattering contributions (direct and diffuse) associated with single and multiple scattering collisions, respectively. In that regard, ICA deals with the partition of a mixed-signal that is composed of linear contributions terms. Unlike principal components, ICA is assuming latent variables having 'non-Gaussian' distributions and independent components that are not necessarily orthogonal. Another important assumption of ICA is the independence between parameters that originated the mixed signal. The initial step of ICA is the whitening (i.e., the covariance of decorrelated variables is an identity matrix) of the original data [15]. The calculation of ICA signal sources is performed by rotating the whitened matrix to minimize the Gaussian distribution behavior ('Gaussianity') between variables. Indeed, the central limit theorem states that a mixedsignal is expected to be more Gaussian than its members.
In this study, ICA calculations were performed by assuming two receivers (i.e., onand off-axis) for measuring two mixed LiDAR signals (X mix 1 and X mix 2 ) with a variable contribution of direct and diffuse backscattering components: where S 1 and S 2 are original signal sources with a dominant direct and diffuse backscattering contribution, respectively, a 1,1 , a 1,2 , a 2,1 , and a 2,2 are weighting coefficients that are influenced by the relative position and orientation of the receivers, and the type of photon interactions with the optical medium. The ICA analysis was restricted to three specific times (hereafter time bin = 110, 160, and 270 or time lags of 55, 80 and 135 ns with respect to the receiver, respectively) representing the leading, descending (i.e., exponential attenuation of power), and trailing portions of each waveform ( Figure 3, Table 2). These portions correspond to different 'energy attenuation regions' where a and b b have different contributions to K sys (e.g., b b >> a and a >> b b in leading and trailing portions, respectively). In other words, LiDAR waveforms are sensitive to different IOPs and energy attenuation processes (e.g., the dominance of direct and diffuse scattering in the leading and trailing portions, respectively) depending on range. The importance of scattering components (i.e., direct/diffuse) and IOPs (i.e., a and b b ) for determining K sys is more balanced as X mix exponentially decreases with distance from the receiver. used to obtain the ICA signal sources (4). The ICA analysis is unable to extract the amplitude of the signal sources, thus the arithmetic mean (μ) and standard deviation (σ) of X1 mix and X2 mix were used to reconstruct the magnitude of S1 and S2, respectively, by applying a z-scores transformation: where IC is the independent component and Sk,λ rec is the reconstructed signal for the signal source S with a dominant scattering contribution k. The proportion of the variability of Sk,λ rec explained by X1 mix and X2 mix was quantified using the coefficient of determination (r 2 ).   Two case studies are presented showing the dominance of 'soft' backscattering features (i.e., n p approximates to n w ) at all detection times ( Figure 3a) and the signal perturbation due to 'hard' backscattering features (n p >> n w ) at detection times between 100 and 150 ns ( Figure 3b). The ICA algorithm used here, also known as FastICA [15], has four major processing steps: data centering by subtracting the mean (1), whitening of centered data based on singular vector decomposition (2), maximization of 'non-Gaussianity' of whitened mixed signals based on kurtosis (3), and normalization/decorrelation of weights used to obtain the ICA signal sources (4). The ICA analysis is unable to extract the amplitude of the signal sources, thus the arithmetic mean (µ) and standard deviation (σ) of X mix 1 and X mix 2 were used to reconstruct the magnitude of S 1 and S 2 , respectively, by applying a z-scores transformation: where IC is the independent component and S rec k,λ is the reconstructed signal for the signal source S with a dominant scattering contribution k. The proportion of the variability of S rec k,λ explained by X mix 1 and X mix 2 was quantified using the coefficient of determination (r 2 ).

Relationships between K sys and IOPs
Optical measurements derived from the ac-9 sensor were useful to interpret the scattering processes determining the LiDAR attenuation coefficient and the identity of LiDAR-derived scatterers. Indeed, the magnitude of K sys varies between c and K d or the diffuse attenuation coefficient of downwelling irradiance [4,18,19]. As c and/or FOV decreases, the light field is expected to be less diffuse and K sys tends to be c. Conversely, K sys tends to be K d when the water becomes more turbid, the FOV becomes larger, and/or the laser beam divergence increases. The final outcome of reducing K sys to K d is a greater contribution of multiple photon collisions to total scattering and a light field that is more diffuse [19]. For each wavelength and detection geometry, K sys values were computed for on-and off-axis receivers (hereafter K sys (on-axis) and K sys (off-axis), respectively) as the slope of log-transform (e-base) backscattering power (X mix ) as a function of range (z): where A L is a constant related to the LiDAR system and n w is the mean refractive index of seawater (i.e., 1.44). The range used for K sys calculations was always within the exponential decay phase of the waveforms and differed between channels (i.e., 8.5-10.7 m and 10.03-12.3 m for onand off-axis, respectively). For each K sys estimate, the mean X mix was computed based on the arithmetic average of the first 100 waveforms (i.e., 1 capture). The slope of each averaged waveform subset was derived by applying a linear regression model type-I [20] to X mix changes as a function of range in m. Waveforms with the presence of strong backscattering features disrupting the exponential decrease of the LiDAR power with range were excluded from the regression analysis. The comparison of LiDAR attenuation coefficients with IOPs-derived from ac-9 measurements was made based on the arithmetic mean of a and c determinations within the water depths matching the LiDAR vertical range used for estimating K sys . To avoid negative backscattering values, an offset of +100 was added to all signals (on-and off-axis) and the signal-to-noise ratio (S/N) was computed as the ratio of X mix (on-axis)/X mix (off-axis). The response of K sys due to changes on IOPs was quantified by the coefficient of determination (r 2 ) after adjusting a linear regression model type I.

SDI
The 2-D structure of Oculus-derived backscattering provides unique information regarding temporal changes on scatterers' distributions that can be mainly attributed to the re-location of optical features due to passive or active motion. Notice that these changes may be associated with variations in backscattering intensity and/or blue/green ratios. The Structural Similarity index (SSI) is a technique widely used in image processing for measuring the similarity between images or 2-D matrices [13]. In our case, the 2-D array is a lidargram or matrix composed of 100 consecutive LiDAR waveforms. Thus, it was assumed that local temporal variability was small with respect to spatial changes of X mix along the boat sampling track and as a function of water depth. Thus, the magnitude of SSI is representative of 200 waveforms (i.e., 100 profiles per lidargram) or 2 s (~5.1 m along the boat direction) and is computed for each element i,j of the lidargram (i.e., horizontal and vertical component, respectively). For each i,j element corresponds to the anomaly of each waveform computed with respect to the median of backscattering values between time bin 110 and 250 (i.e., range = 6.7-14 m). The SSI index was computed as the product of three metrics (luminance,lum, contrast,ct, and structure,st) that apply to each element i,j of lidargrams to be compared (i.e., L1 and L2): where k is the capture time during the survey, α, β and γ are weights set to 1, σ 2 is the variance of element x (i.e., i or j) and σ i,j 2 is the covariance between element i and j, respectively. c 1 , c 2 , and c 3 are constants used to avoid a very small denominator and are affected by the dynamic range. SSI is affected by the size and type of the local window used to smooth the lidargram. In our study, the dynamic range was 256 and the local window was a Gaussian low-pass filter with a size of 11 and a standard deviation of 1.5.
The absence or presence of scattering features between lidargrams was quantified based on two parameters derived from SSI: where n i and n k are the sum of i and j elements, respectively, <SDI> is the arithmetic average of the structural dissimilarity index. For each wavelength, SDI ct in Equation (11) is equivalent to the relative contrast of the backscattering signal between on-axis and off-axis receivers. The range of values for <SDI> and SDI ct is 0-1.

The Spectral Slope of LiDAR Backscattering
The benefits of using a LiDAR with multiple wavelengths in oceanographic applications have been already discussed by Gray et al. [21]. Given their spectral emission and wavelength-dependency on water optical composition, the penetration depth of these systems can be optimized in environments with variable turbidity. Likewise, LiDARs having a spectral resolution allows the identification of scatterers in terms of second-order properties (e.g., the mineral content of particulates). This later advantage was explored here by calculating the spectral slope of LiDAR backscattering (m k ): where k is the receiver (on-or off-axis) with a centered wavelength λ. Notice that m k varies with range or time, thus spectral slopes were analyzed at those time bins described in Section 2.3 and encompassing different 'energy attenuation regions' along the waveforms. Expression (12) was derived by applying a log-transformation to a hyperbolic function proposed for modeling IOPs [22].

Scattering Processes and Shape of Waveforms
Examples of waveforms obtained by Oculus at two wavelengths were shown in Figure 3. The common backscattering volume peak associated with the 'blue' and 'green' on-axis channels was not totally coincident (i.e., photons arriving sooner at the 'blue' receiver). This shift was attributed to the asymmetry of viewing sensor angles making photons arrive first at off-axis receivers (25.2 and 27.3 ns for 'blue' and 'green' channels, respectively) (Figure 3a). Notice that this arrival time difference varied with the water optical properties (e.g., turbidity) and the presence of 'large-sized' scatterers (i.e., 'strong' backscattering features present in on-and off-axis at the leading portion of the waveforms (e.g., 20.3 and 22.4 ns, respectively, Figure 3b). Despite the existence of these 'high' scattering events and differences between detectors in terms of gain and dynamic range (e.g., larger and wider for the green channel), no sensor saturation effects were observed. In a logarithmic space, the beginning of the leading portion of waveforms for on-and off-axis measurements was commonly visualized at 50 ns. For on-axis signals, the exponential decay phase was extended up to 210 ns after which the tail was characterized by a change of slope due likely to a greater contribution of multiple scattering. Conversely, time-resolved backscattering signals for off-axis waveforms were extinct (i.e., S/N < 1) earlier (~195 ns). Perturbations on LiDAR backscattering measurements by a 'hard' scatterers can be seen as a large bulge on the arriving signal (see the second peak at time bin 110 in Figure 3b). In general, the overall impact of this disturbance was related to an increase of signal attenuation after the scattering event and subsequent backscattering oscillations at longer arrival times due to larger contributions of photons having multiple collisions. These oscillations must not be confused with the 'ringing' of the signal (i.e., periodic noise variations after a backscattering saturation event).
The discrimination of time-resolved diffuse and direct backscattering photon contributions to X mix was investigated here based on a subset of waveforms obtained during different times of the day and distances to the shore. As expected, the raw signal for on-axis and off-axis channels was associated to direct-and diffuse-dominated backscattering contributions, respectively (Figure 4). In general, the signal reconstruction was larger (i.e., larger Remote Sens. 2021, 13, 2475 9 of 18 explained variability by the response variable) at time bin 160 followed by time bin 110 and 270 (r 2 up to 0.99, 0.98, and 0.03 with p < 0.001). In the leading portion of the waveforms, the signal reconstruction of direct backscattering returns was higher with respect to that associated to diffuse photon contributions (e.g., explained variability difference up to 13% at λ = 473 nm, Figure 4a,b), and this difference decreased at longer λ. At intermediate detection times (i.e., time bin 160), the direct and diffuse backscattering contributions to S rec were comparable and not influenced by spectral changes. The ICA reconstruction of direct and diffuse backscattering components was at the tail of the waveforms was very poor or null (Figure 4e,f).  The spatio-temporal variability of ICA components is depicted in Figure 5. In general, ICA values for direct and diffuse scattering contributions were less variable during morning hours (i.e., shots 1-800). Also, ICA suggested that 'direct and diffuse photons' covaried positively in the leading portion of the waveforms (Figure 5a,b), a phenomenon that was The spatio-temporal variability of ICA components is depicted in Figure 5. In general, ICA values for direct and diffuse scattering contributions were less variable during morning hours (i.e., shots 1-800). Also, ICA suggested that 'direct and diffuse photons' covaried positively in the leading portion of the waveforms (Figure 5a,b), a phenomenon that was no longer observed at larger distances from the receiver.

Response of Ksys to IOPs
The attenuation of the LiDAR backscattering signal as a function of the range of onaxis waveforms was substantially influenced by changes in water optical properties (Figure 6). This influence was more remarkable with c and within the blue spectral range (p =  = 473 nm (left panels, a,c,e), λ = 532 nm (right  panels, b,d,f), time bin 110 (top panels, a,b), 160 (middle panels, c,d) and 270 (bottom panels, e,f).

Response of K sys to IOPs
The attenuation of the LiDAR backscattering signal as a function of the range of on-axis waveforms was substantially influenced by changes in water optical properties ( Figure 6). This influence was more remarkable with c and within the blue spectral range (p = 0.57, one-tailed t-Student, p = 0.004, Figure 6a). At λ = 473 nm, the a coefficient only explained one-third of attenuation changes on LiDAR backscattering (r 2 = 0.30, p = 0.005, Figure 6b) and no clear relationships were established at the longer wavelength (r 2 = 0.07, p = 0.126). Statistical relationships between K sys (off-axis) and a values were weaker with respect to c-K sys comparisons made at λ= 473 nm and 532 nm (p > 0.05). Covariations between K sys (on-axis) and K sys (off-axis) values were present (r 2 up to 0.32) for waveforms measured within the 'blue' and 'green' spectral range (p = 0.004 and 0.003, respectively, Figure 6c). In general, K sys (on-axis) was larger than K sys (off-axis) (twice on average and up to 2.7-fold at λ = 473 nm).
Morning LiDAR measurements were performed in waters having relatively low turbidity as inferred from c values (e.g., c(488) and c(532) up to 0.34 m −1 and 0.32 m −1 , respectively). Conversely, during the noon and afternoon surveys (i.e., those in shallower areas), the water turbidity was higher (i.e., c(488) and c(532) up to 0.38 m −1 and 0.37 m −1 , respectively) and lidargrams between consecutive captures were less alike. As expected, the mean structural dissimilarity values of waveforms obtained by off-axis receivers were relatively low (<SDI> range = 0.05 ± 0.03, arithmetic mean ±2 standard errors) with respect to those computed for on-axis measurements (t-Student = 13.46 and 11.55 for λ = 473 nm and 532 nm, respectively, two-tailed, p < 0.001, N = 22, Figure 7a,b). Similar to on-axis receivers, off-axis measurements did not show clear spectral differences on <SDI> as values for each wavelength were within the lower and upper bounds of two standard errors (i.e., 95% confidence level). Consistent with <SDI> variations, the SDI contrast between off-axis and on-axis signals tended to have relatively high and more variable values during the noon-afternoon hours (i.e., as high as 0.35) even though these changes were not statistically significant (e.g., SDI ct (morning) vs SDI ct (afternoon), t-Student = −1.82, two-tailed p = 0.083) (Figure 7c,d).

Spectral Slopes of LiDAR Backscattering
The probability distribution function (PDF) of m k values for on-and off-axis measurements obtained during the whole survey and at different detection times are presented in Figure 8. In general, the magnitude of the spectral backscattering slope for on-axis waveforms was substantially larger (i.e., more 100-fold for some samples) and more variable with respect to those corresponding to off-axis measurements.

Spectral Slopes of LiDAR Backscattering
The probability distribution function (PDF) of mk values for on-and off-axis measurements obtained during the whole survey and at different detection times are presented in Figure 8. In general, the magnitude of the spectral backscattering slope for on-axis waveforms was substantially larger (i.e., more 100-fold for some samples) and more variable with respect to those corresponding to off-axis measurements. Most mk values with a probability higher than 50% varied from −4.7 to +6.4 and from −0.5 to +1.3 for on-and off-axis receivers, respectively. For on-axis measurements, the distribution of mk values approximated a Gaussian function for time-bin 110 and 270 ( Figure  8a). Conversely, the normalized PDF of spectral backscattering slopes for time-bin 160 was clearly left-skewed. Likewise, the magnitude of mk values at time bin 160 indicated a decrease of backscattering at longer wavelengths. For off-axis measurements, the shape of the normalized PDF was comparable between different detection times (Figure 8b). For time bins 110, 160, and 270, the arithmetic averages of mk during afternoon surveys were Most m k values with a probability higher than 50% varied from −4.7 to +6.4 and from −0.5 to +1.3 for on-and off-axis receivers, respectively. For on-axis measurements, the distribution of m k values approximated a Gaussian function for time-bin 110 and 270 (Figure 8a). Conversely, the normalized PDF of spectral backscattering slopes for time-bin 160 was clearly left-skewed. Likewise, the magnitude of m k values at time bin 160 indicated a decrease of backscattering at longer wavelengths. For off-axis measurements, the shape of the normalized PDF was comparable between different detection times (Figure 8b). For time bins 110, 160, and 270, the arithmetic averages of m k during afternoon surveys were more positive (i.e., weaker decay of X mix with range) with respect to those computed during morning hours (Table 3). In general, the S/N values for m k calculations were higher at time bin 160 followed in descending order by time bin 110 and 270.

Discussion
The interpretation of results is organized in four main sections encompassing the following topics: the advantages of using ICA for estimating direct and diffuse LiDARderived scattering components (1), the physical meaning of K sys in terms of IOPs (2), the impact of environmental conditions (i.e., turbidity and water depth) on temporal and spatial variability of different types of LiDAR-derived scatterers (3), and interpretation of spectral backscattering variations in terms of size distribution and composition of suspended particulates based on published studies (4).

Direct/Diffuse Backscattering Components
The ICA algorithm was a useful technique to separate direct and diffuse scattering contributions at different time bins and wavelengths. Overall, ICA is a faster and more accurate technique for quantifying scattering sources than traditional Monte Carlo (MC) [9] and PCA (Principal Component Analysis) [23], respectively. Indeed, ICA does not need to follow the trajectory of each photon to elucidate its origin as MC simulations do and, unlike PCA methods, are capable to uncouple correlated interactions by discriminating different probability distributions as derived from higher moments around the mean. In this study, direct and diffuse backscattered photons were mainly associated with signals detected by on-and off-axis receivers, respectively. This can be explained by the larger proportion of photon collisions inside and outside the FOV, respectively. In terms of detection times, the largest reconstruction of direct and diffuse backscattering components corresponded to the exponential decay portion of waveforms where the S/N of backscattered photons was higher with respect to those values characteristics of the leading and tail waveform sections. Also, the better discrimination of direct and diffuse components at time bin 160 was attributed to the greater 'Gaussianity' of probability distribution functions for photons arriving at time bin 110 and 270 (i.e., where multiple scattering contributions increase).

LiDAR vs Ac-9 Optical Properties
The exponential decay of LiDAR backscattering power with range showed systematic differences among waveforms captured by off-axis and on-axis receivers. In general, the signal attenuation was more remarkable for on-axis waveforms and this was attributed to the dominance of single scattering (backscattering + forward-scattering). As multiple scattering increases and the light field becomes more diffuse (i.e., off-axis measurements), the magnitude of K sys decreases approaching K d [4,18,19]. For on-axis measurements, K sys had a stronger covariation with c values and suggests that scattering rather than absorption is the driving process modulating K sys in our measurements. As expected, this effect was more pronounced in the blue spectral range where backscattering of particulates and water generally increases with respect to green wavelengths.
The apparent insensitivity of K sys (off-axis) to changes on IOPs distributions was likely attributed to the uncoupling of two important factors determining multiple scattering: path-length and water turbidity. For a constant c, multiple photon collisions are anticipated to augment at longer distances from the receiver. This multiple scattering regime is usually identified as a change on the backscattering slope (i.e., ln(X mix ) as a function of range). As c increases, this slope break is expected to move with the maximum path-radiance or common backscattering volume toward earlier detection times [23]. In fact, additional backscattering due to higher c values and greater contribution of multiple scattering leads to a shape deformation in the original waveform that is mainly characterized by backscattering slope changes in the leading and trailing sections and concurrent non-linear variations on K sys . The optical configuration of the Oculus system is consistent with K sys values that approximate c rather than K d . Indeed, the composite product cR [19], where R is the illuminated spot at a specific range (e.g., 6.2, 9, and 15.1 m for time bins of 110, 160, and 270, respectively) was always very small (0.07 to 0.19) compared to 1. Likewise, light scattering was the dominant process determining LiDAR backscattering attenuation in our study and explained two important findings: (1) the response of K sys (on-axis) to c changes was larger with respect to that associated with a variations, and (2) the larger magnitude of K sys (on-axis) with respect to K sys (off-axis). In the last case, the difference between K sys values suggests that K sys (on-axis) has an additional attenuation term due to scattering as K sys (off-axis) is mainly driven by light absorption.

Diel and Spatial Patterns of Scatterers
A consistent pattern revealed by ICA at all detection times was the relatively low variability of ICA components during morning hours. This phenomenon was likely related to the sampling design and environmental differences regarding water types. In the first case, morning surveys were part of a time series and explain why transects (i.e., noonafternoon profiles) were characterized by having larger changes on ICA values as spatial measurements include two sources of physical variability (local + advective). In the second case, the water optical properties of the morning dataset were different from those measured during the noon-afternoon datasets. Indeed, the water turbidity as inferred from c suggested a predominant oceanic (coastal) water type during morning (noon-afternoon) profiles. Since the average size of particulates increases with turbidity [17,24], larger backscattering features were likely more abundant late during the day when surveys were closer to the shore.
Assuming a negligible spatial variability between captures, the comparison of <SDI> and SDI ct values for different sampling locations and times of the day suggested changes in mobility of scatterers during our surveys. Indeed, maximum values of similarity between lidargrams were found offshore. This observation was confirmed based on contrast index calculations (i.e., maximum SDI ct values near the coast). The lack of coherence between LiDAR backscattering profiles near the coast was associated with higher water turbidity levels as inferred from c and associated changes in particle characteristics as discussed above. However, the motion of large-sized scatterers due to active swimmers (e.g., fish) [25] was likely another influencing factor. Reef fish along the south Florida shelf is known to be highly aggregated near the coast [26], thus it is likely that observed LiDAR backscattering patterns were partially related to fish distribution differences across the shelf. In general, <SDI> values associated with off-axis measurements were smaller with respect to those derived from on-axis measurements. This is not unexpected as off-axis measurements are dominated by diffuse scattering contributions that are less influenced by the presence of relatively 'strong', 'large-sized', and less common backscattering features (e.g., jellyfish) [27]. These relatively rare optical features were in part responsible for augmenting <SDI> associated with on-axis waveforms and preferentially those obtained nearshore and at shorter wavelengths.

Spectral Backscattering Variations
In general, on-axis waveforms were characterized by having a greater proportion of positive m k values (i.e., X mix increases at longer wavelengths) with respect to those derived from off-axis measurements. This phenomenon was attributed to the greater sensitivity of on-axis measurements to 'large-sized' optical features (i.e., geometric cross-section much larger than LiDAR wavelength). Time-series off the Massachusetts coast have shown consistent differences in spectral backscattering between different particle size distributions as derived from backscattering meters [14]. Indeed, Slade et al. found more negative spectral backscattering slopes (λ = 488-715 nm) when suspended particulates within the size range 5-50 µm were dominated by finer size fractions. Loisel et al. [28] computed the spectral backscattering slopes for different marine regions around the globe based on satellite observations and concluded that most estimates vary between 0 and −3.5 with more positive values associated with eutrophic zones where suspended particulates are larger. Lastly, tank experiments using LiDAR [21] found that spectral backscattering slopes of Arizona dust (mean diameter = Dm = 4.5 µm) were generally smaller (−1.72 to 0.57) with respect to those (−0.24 to 2.17) derived from large-sized organic particles associated to a phytoplankton culture of I. galbana (Dm = 6.5 µm). Notice that LiDAR measurements made in [21] correspond to a biaxial geometry (i.e., detector and source are non-collocated), thus their resulting waveforms were more alike to our off-axis backscattering determinations.
The spectral composition of backscattered photons differed between relatively shallow and deep water locations (i.e., predominantly coastal and oceanic conditions, respectively) as m k values become more positive closer to the shore and during the afternoon (Table 3). This spatial trend was likely related to the greater proportion of 'large-sized' particulates and higher water turbidity near the shore. Spatial patterns on spectral beam attenuation, an optical proxy for particle size distribution [14], supports this hypothesis as high c(488)/c(532) ratios (i.e., a greater contribution of 'small-sized' particulates) tended to decrease closer to the shore (i.e., morning samples) ( Figure A1, Appendix A). Gray et al. [21] pointed out a substantial increase of m k (e.g., −1.7 to −0.24 and 0.57 to 2.17 in clear and turbid waters, respectively) in turbid waters (i.e., up to a 6.5-fold increase of c(550)). Also, results on [21] showed composition effects on m k with more positive and negative values associated with organic-dominated (−0.24 to +2.17) and mineral-dominated (−1.72 to +0.57) particulates, respectively. In this study, the range of m k values derived from off-axis receivers suggests that particle assemblages had an intermediate chemical composition between inorganic-rich and organic-rich case studies. Likewise, the increase of m k values toward the coast in our surveys indicates that particle composition effects on m k were secondary with respect to those associated with PSD and/or turbidity changes.

Conclusions
The discrimination of underwater optical features and characterization of scatterers by standard LiDAR configurations is limited due to their relatively poor spectral resolution and low signal/noise ratios. In this study, Oculus, a new multispectral LiDAR system was applied to understand scattering sources and their relationships with IOPs and scatterer types in shelf waters off the South Florida coast. In general, ICA suggested a more defined separation between direct and diffuse scattering contributions along the exponential decay of X mix due to a higher S/N and the 'non-Gaussianity' behavior of the probability distribution function. This waveform region is commonly used to compute K sys (range or non-range resolved LiDAR attenuation coefficient) and estimate IOPs. In our case, K sys variability was linked to c as supported by in situ measurements and theoretical considerations. The complimentary use of SDI/SDIct and m k was useful to identify scatterers in terms of their properties and distribution patterns. Indeed, structure dissimilarity indexes suggested a greater mobility of scatterers near the coast where m k also indicated a greater predominance of relatively large-sized particulates and more turbid conditions.