Automated Riverbed Sediment Classification Using Low-Cost Sidescan Sonar

Classifying subaqueous riverbed sediments by grain size provides a means to parameterize important boundary conditions for studies of channel hydraulics, aquatic ecology, sedimentation, geomorphic change, and sediment transport. Sampling riverbed sediment deposits with sufficient spatial density and coverage is particularly difficult when the water is too turbid or deep to image the bed from aerial platforms, or too swift or deep to wade to obtain physical samples. Conventional underwater imaging systems are limited by light attenuation, turbidity, and the small spatial footprint of images (Buscombe et al. 2014a). Sidescan sonar has the potential to address this technical shortfall (Kaeser and Litts 2010) by recording acoustic signals that can be analyzed in a manner that distinguishes varied substrate textures on the basis of related patterns of echo intensities.

Sidescan sonar, deployed from a moving vessel, produces photograph-like images of river and bed texture (Fig. 1). The transducer sends out a high-frequency (typically several hundred kHz) acoustic beam perpendicular to the vessel heading on either side (port and starboard) and records the amplitude of the returning echoes from a wide swath (Blondel 2009). One ping constitutes the simultaneous acquisition of data from the two sidescan beams at an instant, returning a swath composed of pixels whose intensity relates to the echo strength, determined by acoustic impedance and reflection at those locations. A small strip of the bed is imaged with each ping, building an echogram that provides near continuous coverage as the vessel moves along-track (up or downstream). Spatial distributions of bottom texture visible within echograms are related to the roughness of the bed, which provides a means to map bed-sediment types (Collier and Brown 2005). Sidescan pixels are typically of order centimeter to decimeter in size, therefore do not resolve individual sand and gravel clasts but do resolve cobbles and boulders, yet the texture in the sonar echogram for patches of sand and gravel should be distinct and quantifiable. In this paper, texture is defined as the frequency of change in arrangement of small-scale surface roughness, which in turn is a measure of the statistical variation in the distribution of bed sediment grain sizes. Texture is used as a 2D quantity without formal definition but is, in concept, closely related to autocorrelation. Less textured surfaces (i.e., over smooth beds) in sonar echograms show less variation in adjacent sidescan pixels over space; they are less textured and have smaller texture lengthscales. In contrast, rough beds in sonar echograms have more texture because they have relatively large spatial variations in adjacent pixels over space (caused by reflections and shadows): the spacing between textural elements is larger and the rough bed is characterized by larger average texture lengthscales.

Several inexpensive sidescan units that are easily mounted to small vessels are commercially available and suitable for substrate imaging. These low-profile systems have a minimal draft requirement so they are especially suitable for imaging in very shallow water. Lightweight and with low power demands, they can be operated by one person in any river or stream navigable by a small boat. Kaeser and Litts (2010) reported that one such system (the sidescan sonar on the Humminbird fishfinder) was of sufficient quality for bed imaging in shallow rocky streams. Kaeser et al. (2012) produced bed-sediment maps by merging and superimposing imagery onto a base map in a Geographical Information System (GIS), identifying regions of similar bed texture by visual inspection. In contrast, multibeam sonar systems are typically heavier, more expensive, with greater computational and power demands, and require a greater draft as well as specialist knowledge to operate.

Visual interpretations of bed texture are labor intensive and therefore impractical for mapping long river segments, whereas automated sediment classification from echograms would be an objective and rapid means to do so. Classification of bed sediments by characterizing the spatial variability of echogram textures using frequency-domain analysis (Ruet et al. 1985; Pace and Gao 1988; Tamsett 1993; Atallah et al. 2002) has a number of advantages: it is computationally efficient; makes no assumptions about the particular sidescan device used or the nature of scattering of sound by various substrates; and does not require an acoustically calibrated system (Ruet et al. 1985).

This study suggests a method that uses wavelet spectral techniques for bed sediment classification. The approach relies on computation of average textural lengthscales based on the variance associated with different wavelengths of bed roughness elements and their acoustic shadows. The lengthscales may be used in an objective of classification of sediments based on their qualitative Wentworth class (silt, sand, gravel, cobble, and boulder), commonly used in studies of aquatic ecology (Cross et al. 2013), sediment budgets (Grams et al. 2013), and models for sediment transport (Wilcock and Kenworthy 2002).

Unlike Fourier analyses, the requirement of signal stationarity is not so strict for wavelet-based analyses (Torrence and Compo 1998), making it more suitable for short, intermittent spatial series with non-normal distributions (Baas 2006) typical of textures of natural materials (Qi and Neupauer 2008; Buscombe 2013). The wavelet transform has long been used to classify texture in images (Charalampidis and Kasparis 2002) including sidescan textures (Little and Smith 1996) for these reasons. The authors’ approach facilitates classification at high-resolution by enabling analysis of small nonstationary series of echo levels. Echogram textures are characterized by variance associated with different texture wavelengths and because wavelets are well-resolved in both the spatial and frequency domains, in contrast to an inherent trade-off between the two in Fourier-based spectral analyses (Torrence and Compo 1998), averaging over space provides so-called global spectral estimates (Farge 1992).

The data recorded by a Humminbird 998c unit include port and starboard sidescan echograms (at 455 kHz), and echograms from two downward-looking echosounders (at 83 and 200 kHz, respectively) measuring nadir depth, position (Global Mercator coordinates), and vessel speed. Sound celerity calculations and range-corrections are carried out in real time using measured water temperature assuming zero salinity and homogeneous mixing.

The system was evaluated through deployments on the Colorado River in Grand Canyon National Park, Arizona, in May 2012 and May 2014, and the Penobscot River in Maine during September 2013. On the Colorado River, the system was used to image pools between rapids and riffles, throughout a 50 km segment of the river. The sonar was mounted to a 4.57 m long motorized inflatable cataraft away from sources of turbulence such as the free surface, motor, and strakes. In both settings, the draft to transducer face was 20 cm, and vessels were allowed to drift slowly downstream during data collection at a nearly constant speed and heading. Therefore, corrections for instabilities in vessel attitude were unnecessary. In Maine, two reaches were imaged from a 4.26 m long inflatable raft, one between rapids on the West Fork of the Upper Penobscot River above the town of Millinocket and one immediately downstream from the head of tide on the river near the city of Bangor.

Both settings present challenges to the sensing of bottom conditions. The Colorado River can be characterized by relatively high suspended concentrations that can attenuate light and high-frequency sound; however, the concentrations were low ($<100\mathrm{mg}/\mathrm{l}$) during the time of the survey. The Penobscot has characteristically low suspended sediment concentrations ($\ll 100\mathrm{mg}/\mathrm{l}$), and visualization of bottom conditions is impeded by relatively high concentrations of natural tannins. Both riverbeds are dominated by a mix of sand and gravel with insignificant deposits of cohesive sediments.

The Colorado River bed in Grand Canyon is highly dynamic, encompassing the full range of noncohesive sediments, deposited in complicated spatial patterns (Grams et al. 2013). Tributary inputs of sand and gravel are variable, and hydraulics of individual pools and eddies vary considerably in space and in response to varying dam operations. The bed of the Penobscot (Fig. 2) is generally more uniform, predominantly embedded gravels with small isolated sand patches in sections with modest or low gradients (Dudley and Giffen 2001). Isolated patches of cobbles, boulders, and bedrock are present in the lower reaches of the river near locations of recent dam removal projects but of more limited extent below the head of tide on the river (PRRPS 2008). Aggregations of coarse materials often correspond to locations of abrupt bed elevation drops in the Upper Penobscot. Rapids were not imaged as part of the surveys of either river, owing to the deleterious effects of excessive aeration on the high-frequency acoustics.

A generalized routine is described for estimating the average lengthscales of acoustic fluctuations in signals from sidescan sonar, which are the alternating patterns of high and low backscattered sound from objects on the bed and the shadows they cast. The lengthscales of fluctuation should scale with the average size of textural elements on the bed. The average lengthscale of signal fluctuation is calculated using a Morlet wavelet transform, and then the bias caused by varying sonar geometry is removed.

Let $I$ be a 2D portion of echogram collected by the sidescan as the boat moved downstream [Figs. 3(a and e)], composed of across-track (river spanwise) scans composed of pixels distributed on the bed. The across-track scans are concatenated in the along-track (streamwise) direction, such that $I$ is composed of $N$ rows of discrete scans $y_n$ ($n=0$. $N-1$ sampled at equal spacing, i.e., each pixel). For each $y_n$ [Figs. 3(b and f)], the linear trend is removed, then zero-padded to the next power of 2 greater than $N$. The power autospectral density of each $y_n$ is estimated using a wavelet transform, by convolving each $y_n$ with a scaled and translated version of a normalized mother wavelet function ${\psi }_0$ (Kumar and Foufoula-Georgiou 1997)where $\delta y$ is the spacing in $y_n$; the apostrophe indicates transpose; and the asterisk indicates complex conjugate. By varying the scale, $s$ (wavelength in pixels), and translating in location, $n$ (position in pixels along the sequence), a one-dimensional series is expressed in a two-dimensional parameter space ($n$, $s$) [Figs. 3(c and g)]. Here, the complex-valued Morlet is used for ${\psi }_0$ because it is well resolved in both spatial and frequency domains (Torrence and Compo 1998, Buscombe 2013). The variance in autospectral power across $N$ spectral densities $W_n(s)$ provides a measure of global autospectral variance, as a function of scale (s), computed over multiple locations ($n$) [Figs. 3(d and h)] as

Highly textured echograms [Fig. 3(e)] vary more strongly over multiple scales [Fig. 3(f)], therefore the energy in the continuous wavelet transform is broad-banded [Fig. 3(g)] and the width of the variance spectrum—qualitatively, the number of important wavelengths in the data—increases [Fig. 3(h)]. For echograms with little texture [Fig. 3(a)], the variance spectrum is narrow [Fig. 3(d)] because fewer wavelengths are required to describe the data. The mean textural lengthscale, $w$, computed from the normalized variance spectrum (${\sigma }^2/\sum {\sigma }^2$) using spacing $\delta s$describes the average lengthscale of signal fluctuations caused by backscattering from objects on the bed and their acoustic shadows.

Calculations are more tractable in 1D because an echogram has an irregular spatial grid. Time-varying current speed, vessel speed, and sonar ping-rate (due to variations in water depth and sound celerity) mean that physical distances between successive pings (along-track spacing) are rarely regular. Likewise, temperature-driven variations in sound celerity can cause variations in across-track spacing. Therefore, to convert $w$ from units of pixels to meters in the across-track direction, $w$ is normalized by $\pi /2{\sin }^{-1}(c/tf)$. This scaling, where $c$ is the speed of sound in water, $t$ is the physical length of the transducer array, and $f$ is the frequency of sound, is related to the sonar geometry and beam spreading. In the along-track direction, $w$ is integrated with respect to $\delta x$, the along-track spacing (meters) between successive pings. The echogram is segmented into windows, $I$, each of size $N\times N$ and overlapping by $N-1\times \mathrm{N}-1$ pixels. Texture lengthscale is calculated for each $1\times N$ along-track and across-track series, and results averaged. By assigning the computed value of $w$ to the coordinate at the center of each window, it is possible to build an output matrix that is identical in size (and resolution) to the input echogram. The windowing procedure is essentially the same as that described, in more detail, by Buscombe et al (2014b) for spatially explicit analysis of multibeam sonar backscatter [for a video clip of an annotated worked example of automated texture lengthscale mapping from a wavelet windowing procedure on sidescan echogram data, see Buscombe (2015)].

This procedure provides a measure of the average lengthscale of signal fluctuations and is related to the average lengthscale of bed texture elements. The final step is to remove the bias associated with the inverse relation between the length of the acoustic shadow, $L$, cast by an object on the bed with height $h$ and the angle between the acoustic beam and the bed (the grazing angle, $\theta$), related by (Blondel 2009)where $H$ = height of the sonar above the bed; $R$ = slant-range to the object. The above says that, for an object of given $h$, its shadow will lengthen as grazing angle, $\sin (\theta )=H/R$, decreases. Therefore, for a given $h$, there is a $\cos (\theta )$ bias in $L$ and, under the reasonable assumption that ${\mathrm{L}}_{\text{object}}/{\mathrm{L}}_{\mathrm{shadow}}\approx$ constant over small (order square meter) window $I$, average texture lengthscales are recomputed as $w=w\cos (\theta )$.

The average lengthscale of bed texture elements is a physical criterion used to delineate bottom areas of unique texture within a sidescan echogram. The textural lengthscale is not a direct measure of grain size. Rather, it is a statistical representation that integrates over many attributes of bed texture, of which grain size is the most important.

The method was applied to all data on the Colorado and Penobscot Rivers. The chosen window size was 200 square pixels, which represents up to a few square meters of bed, depending on $\delta y$ (across-track) and $\delta x$ (along-track).

The automated texture map method was evaluated in two ways. First, texture lengthscales on a heterogeneous (sand through to boulder) patch of the Colorado River bed (Fig. 4) were measured on-screen under zoom, by eye, in the three regions shown, and histograms of these measurements were compiled (Fig. 4, inset). The length of textural elements (the high-intensity portions of the textural patterns in between successive shadows) were measured in the along-track and across-track directions. These manual measurements were compared to the equivalent metric from the automated routine. The example shows that the method satisfactorily estimates texture lengthscales for the three regions, which are interpreted to be boulder, coarse gravel, and fine gravel patches, respectively.

Second, texture lengthscales were computed for several scans taken at two sites in Marble Canyon over larger areas (up to several hundred meters longitudinally, and up to full channel width), merged into continuous texture lengthscale maps, and compared with georeferenced underwater video observations of the bed. More details of the sites and ground-truth data are found in Buscombe et al. (2014c), which used similar data to validate a substrate classification technique based on high-resolution multibeam acoustics. At both sites, the bed is predominantly sandy with isolated patches of rock and gravel. At river mile 30 (48 km downstream of Lees Ferry, Arizona), three scans taken in May 2012 over a $\sim 200\times \sim 60\mathrm{m}$ reach were analyzed (Fig. 5), texture lengthscales were gridded at 1 m and compared to 20 georeferenced observations of bed substrate type (Fig. 6). A similar analysis was conducted using three scans taken in May 2014 at river mile 61 (98 km downstream of Lees Ferry) over a $\sim 450\times \sim 60\mathrm{m}$ reach (Fig. 7), and results were compared to 36 bed sediment observations (Fig. 8). In both cases, close qualitative agreement demonstrate the potential use of the technique for automated substrate mapping. In addition, there is reasonable agreement between these substrate classifications and those based on high-resolution multibeam backscatter at the same sites [Figs. 6 and 4 in Buscombe et al. (2014c) for substrate maps at river mile 30 and 61, respectively]. A more thorough validation of the technique will involve detailed description of the objectives of a specific case study and steps taken to meet acceptable levels of accuracy and precision in substrate classifications.

A spectral method has been proposed to objectively create continuous maps of bed texture lengthscales from sidescan echograms, for the purpose of spatially-explicit characterization of bed sediment. An overlapping windowing procedure is performed, and a direct statistical measure of the average lengthscale of sediment textures on the bed is calculated per window. In situations where a scaling relationship exists between this average texture lengthscale and sediment type, this technique should find utility. Sensitivity to window size, and its optimal size, will be the subject of further study. Aggregating (by averaging, weighted averaging, or otherwise) results from analyses using several window sizes might reduce sensitivity in results through regression toward the mean. It might also be possible to determine optimal window size statistically from the data.

In agreement with previous studies (Kaeser and Litts 2010), the authors observed that the Humminbird sidescan device produces signal returns of sufficient quality for mapping substrates in gravel and sand-gravel riverbeds. The effects of organic materials (e.g., aquatic vegetation and woody debris) and deposited cohesive sediments have not yet been evaluated. Further work is also required to assess the applicability of the method for echograms degraded by highly aerated water or high concentrations of suspended sediment.

The technique proposed in this paper extends previous spectral methods in four ways. First, the classification is spatially explicit, i.e., on a pixel-by-pixel rather than on a swath-by-swath basis (e.g., Atallah et al. 2002) or over larger areas (Ruet et al. 1985; Pace and Gao 1988; Tamsett 1993). In rivers it is often unreasonable (e.g., Nelson et al. 2009) to assume that sediments will not vary across scales of tens to hundreds of meters. Second, the signal-processing of textures directly results in a classification without any intermediary statistical steps such as those based on clustering or machine learning. Third, the method quantifies textures in physical units of length that correspond to surface roughness, rather than acoustic (Collier and Brown 2005) or nondimensional, non-physical quantities (Cochrane and Lafferty 2002; Huvenne et al. 2002; Atallah et al. 2002). Finally, simple consideration is given to sonar geometry, correcting spectrally computed lengthscales for any shadow-lengthening bias introduced by low grazing angles.

Thanks to dedicated field technicians and river guides. The manuscript benefited greatly from comments by the associate editor and two anonymous reviewers. This work was funded by the Glen Canyon Dam Adaptive Management Program administered by the U.S. Bureau of Reclamation. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the United States government. Tools for carrying out the analyses described in this paper have been developed in an open-source, cross-platform language (Python), and made computationally efficient using parallel processing and an optimising static compiler (Cython). This software is freely available under a GNU license at: https://github.com/dbuscombe-usgs/PyHum.