A Preliminary Assessment of a Newly-Defined Multispectral Hue Space for Retrieving River Depth with Optical Imagery and In Situ Calibration Data

Abstract

Bathymetry is a key element in the modeling of river systems for flood mapping, geomorphology, or stream habitat characterization. Standard practices rely on the interpolation of in situ depth measurements obtained with differential GPS or total station surveys, while more advanced techniques involve bathymetric LiDAR or acoustic soundings. However, these high-resolution active techniques are not so easily applied over large areas. Alternative methods using passive optical imagery present an interesting trade-off: they rely on the fact that wavelengths composing solar radiation are not attenuated at the same rates in water. Under certain assumptions, the logarithm of the ratio of radiances in two spectral bands is linearly correlated with depth. In this study, we go beyond these ratio methods in defining a multispectral hue that retains all spectral information. Given n coregistered bands, this spectral invariant lies on the $(n-2)$-sphere embedded in ${\mathbb{R}}^{n-1}$, denoted ${\mathbb{S}}^{n-2}$ and tagged ‘hue hypersphere’. It can be seen as a generalization of the RGB ‘color wheel’ (${\mathbb{S}}^1$) in higher dimensions. We use this mapping to identify a hue-depth relation in a 35 km reach of the Garonne River, using high resolution (0.50 m) airborne imagery in four bands and data from 120 surveyed cross-sections. The distribution of multispectral hue over river pixels is modeled as a mixture of two components: one component represents the distribution of substrate hue, while the other represents the distribution of ‘deep water’ hue; parameters are fitted such that membership probability for the ‘deep’ component correlates with depth.

1. Introduction

Characterizing river channel bathymetry by average parameters (e.g., average depth, cross-section area, hydraulic radius) or sparsely distributed cross-sections may be sufficient for certain applications; however, a precise description is needed for applications in geomorphology or ecology, where local flow parameters are known to vary at the scale of a few times the full bank width in the streamwise direction with alternating morphological units such as pools, runs, and riffles [1,2].

Dense and precise bathymetric survey requires active measurement techniques such as bathymetric LiDAR using water-penetrating wavelengths in the visible spectrum [3,4], or acoustic soundings [5]. Unfortunately, these techniques rely on an artificial illuminating source (either electromagnetic or acoustic) and thus have a high operating cost and are very difficult to use over large areas.

A very interesting alternative is offered by passive optical methods, which have been used in marine, coastal, and more recently river settings since the pioneering work of Lyzenga [6]. Here, natural sunlight is the illuminating source and depth is retrieved by analyzing the radiance signal received at the sensor. According to Marcus and Fonstad [7], this simple method may actually be ‘the only viable method for measuring, monitoring and mapping a large suite of in-channel river parameters continuously at sub-meter resolution’.

Following this introduction, Section 2 of this paper will provide a quick overview of the principles behind remotely-sensed estimation of river bathymetry using spectral ratios between bands in multispectral images. Then, in Section 3, we will temporarily move away from this specific application and propose a formal definition of a more general multispectral hue. Our procedure for estimating depth using this multispectral hue will be the core of Section 4, and results will be presented in Section 5. Then, we will summarize our findings and propose perspectives in Section 6.

2. Principles of Remotely-Sensed Estimation of River Bathymetry Using Passive Optical Measurements

As illustrated in Figure 1, the total upwelling spectral radiance (in W·m⁻²·sr⁻¹) over water in a given spectral band with central wavelength $\lambda$ can be split into four components (e.g., Legleiter et al. [8]):

$$L_{\uparrow }(\lambda )=L_b(\lambda )+L_c(\lambda )+L_s(\lambda )+L_p(\lambda )$$

(1)

• $L_b$ is the radiance originating from the Lambertian (diffuse) reflection from the bed substrate;
• $L_c$ is the volume radiance of the water column of depth h, which is basically sunlight that has been backscattered upwards before reaching the bottom;
• $L_s$ is the surface radiance due to specular reflections at the air–water interface. $L_s$ can make up a large fraction of $L_{\uparrow }$ for certain geometries or viewing angles (‘sun glints’);
• $L_p$ is a path radiance due to atmospheric scattering.

Providing that we can account for both atmospheric path radiance and surface radiance, Philpot [10] derived the following expression for the total radiance measured at a given pixel of the sensor:

$$\begin{array}{ccc}\hfill L_{\uparrow }(\mathbf{x},\lambda )& \approx & L_b(\mathbf{x},\lambda )+L_c(\mathbf{x},\lambda )\hfill \\ & \approx & \underset{R_{\mathrm{mix}}(\mathbf{x},h,\lambda )\hspace{0.17em}‘\mathrm{mixture}’\hspace{0.17em}\mathrm{reflectance}}{\underbrace{\left[\left(R_b(\mathbf{x},\lambda )-R_{c,\infty }(\lambda )\right)e^{-K(\lambda )h(\mathbf{x})}+R_{c,\infty }(\lambda )\right]}}{C}_0{T}_a(\lambda )E_{\downarrow }(\lambda )\hfill \end{array}$$

(2)

The rightmost term $E_{\downarrow }(\lambda )$ is the downwelling solar irradiance in the spectral band (in W·m⁻²), while the product $C_0{T}_a$ (in sr⁻¹) accounts for both atmospheric and air–water interface transmission ($C_0$ can be considered roughly independent of wavelength). These terms are then multiplied by a reflectance term (brackets) which, according to Philpot, can be written as a mixture between two ‘end-members’:

• The first end-member is the reflectance of bed/substrate, $R_b$, which depends on both location $\mathbf{x}$ and wavelength $\lambda$,
• The second end-member is the reflectance $R_{c,\infty }$ that would be measured over an ‘infinitely deep’ water column.

The weight of the bed reflectance $R_b$ in the mixture reflectance $R_{\mathrm{mix}}$ is an exponentially decreasing function of depth, with an attenuation coefficient $K(\lambda )$, which depends on the wavelength (Beer–Lambert law).

While the solar irradiance $E_{\downarrow }(\lambda )$ depends on the time of the day and on cloudiness, taking the ratio of solar irradiances in two different spectral bands cancels these variations. Hence, it is convenient to take the ratio:

$$\frac{L_{\uparrow }(\mathbf{x},{\lambda }_i)}{L_{\uparrow }(\mathbf{x},{\lambda }_j)}=\left[\frac{\left(R_b(\mathbf{x},{\lambda }_i)-R_{c,\infty }({\lambda }_i)\right)e^{-K({\lambda }_i)h(\mathbf{x})}+R_{c,\infty }({\lambda }_i)}{\left(R_b(\mathbf{x},{\lambda }_j)-R_{c,\infty }({\lambda }_j)\right)e^{-K({\lambda }_j)h(\mathbf{x})}+R_{c,\infty }({\lambda }_j)}\right]\frac{C_0{T}_a({\lambda }_i)E_{\downarrow }({\lambda }_i)}{C_0{T}_a({\lambda }_j)E_{\downarrow }({\lambda }_j)}$$

(3)

Two common assumptions are then made to proceed with this expression:

• First, the reflectance $R_{c,\infty }$ of deep water is assumed to be small compared to the bed reflectance $R_b$ (for any substrate and any wavelength). Deep water reflectance is certainly low if the water is clear, but this assumption breaks down in the presence of sediment, dissolved organic matter, algae, etc. Furthermore, the bed might have a very low reflectance. Quartz sand might be highly reflective, but mud and rocks (especially rocks coated with a microbial film) can be quite dark;
• Second, the ratio of bed reflectances at wavelengths ${\lambda }_i$ and ${\lambda }_j$, $\frac{R_b(\mathbf{x},{\lambda }_i)}{R_b(\mathbf{x},{\lambda }_j)}$ is approximated as uniform in space.

With these assumptions, Equation (3) can be simplified and boils down to a linear relationship between depth and the logarithm of the spectral ratio:

$$X_{ij}=ln\frac{L_{\uparrow }(\mathbf{x},{\lambda }_i)}{L_{\uparrow }(\mathbf{x},{\lambda }_j)}\approx \underset{\mathrm{slope}({\lambda }_i,{\lambda }_j)}{\underbrace{\left(K({\lambda }_j)-K({\lambda }_i)\right)}}h+\underset{\mathrm{intercept}({\lambda }_i,{\lambda }_j)}{\underbrace{ln\left(\frac{R_b({\lambda }_i)}{R_b({\lambda }_j)}\right)+ln\frac{C_0{T}_a({\lambda }_i)E_{\downarrow }({\lambda }_i)}{C_0{T}_a({\lambda }_j)E_{\downarrow }({\lambda }_j)}}}$$

(4)

Slope and intercept depend on optical properties of the atmosphere, water, and substrate at each wavelength but, if they can be considered invariant across the image, Equation (4) provides a predictor for depth. It has been used in numerous studies (see Shah et al. [11] for a recent review), either with the calibration of slope and intercept against in situ depth measurement, or without calibration; in the latter case, attenuation coefficients $K(\lambda )$ have to be modeled physically [12,13].

For a multispectral image with n coregistered bands, $n(n-1)$ spectral ratios can be computed, but only $(n-1)$ are independent (keeping those between adjacent bands for instance) since

$$X_{ik}=ln\frac{L_{\uparrow }({\lambda }_i)}{L_{\uparrow }({\lambda }_k)}=ln\left(\frac{L_{\uparrow }({\lambda }_i)}{L_{\uparrow }({\lambda }_j)}\frac{L_{\uparrow }({\lambda }_j)}{L_{\uparrow }({\lambda }_k)}\right)=X_{ij}+X_{jk}$$

(5)

Hence, a ratio–depth relation can be obtained by identifying a single, optimal pair of bands [8], or using several log-transformed ratios given by several pairs [14]. The use of multiple ratios (or other predictors) will increase model flexibility, though attention must be paid to multicollinearity [15]. In fact, retaining all spectral information from the $(n-1)$ independent ratios amounts to defining a multispectral ‘hue’, or ‘pure color’, independent of illumination conditions. The definition of such a spectral invariant is presented in Section 3.

3. Beyond Spectral Ratios: A Definition of Multispectral Hue

3.1. Rationale for the Definition of Multispectral Hue

3.1.1. RGB Case: The Color Wheel

Let us start with an example that will be familiar to many: many computer programs offer a color definition tool that relies on the kind of ‘color wheel’ shown in Figure 2. In this tool, a hue or ‘pure color’ can be selected along a circle, with an angle ranging between 0 and 360°. Two other operations can then be performed to obtain the final red, green and blue components:

• The pure color can be partially desaturated, i.e., mixed with a white component. Such a mixing can occur for example with specular reflections, which will appear nearly white under white illumination if the refractive index of the material is slowly varying with wavelength (this is the case for water, or for an object coated with varnish).
• The overall brightness or value of the color can be decreased. Physically speaking, this other operation is similar to either reducing the intensity of the light source, or changing the orientation of an object’s surface normals so that less light is diffused towards the sensor (‘shaded face’ effect).

3.1.2. Notations

In this section and the following ones, we will use bold symbols to denote points, vectors or matrices in a multidimensional space. Denoting ${\mathbf{v}}^T$ the transpose of $\mathbf{v}$, $[1,1,1]$ denotes a row vector while ${[1,1,1]}^T$ denotes a column vector.

3.1.3. The Reason Why Multispectral Hue Should Have $(n-2)$ Degrees of Freedom with n Bands

Defining the multispectral hue of a pixel rigorously amounts to reversing the two operations described above, that is, resaturating the color and enhancing/normalizing its brightness. Because of these two steps, multispectral hue should have only $(n-2)$ degrees of freedom. This property is obviously verified with $n=3$ bands and the resulting color wheel, which provides only one angular degree of freedom ($n-2=1$).

We start our general definition of multispectral hue with the rationale of Montoliu et al. [16]. Given a pixel whose values in the n available bands form a vector $\mathbf{C}={[C_1,C_2,\dots ,C_n]}^T$ of digital counts, they first define a vector $\mathbf{c}$ with entries

$$c_i=C_i-\min \left\{C_1,C_2,\dots ,C_n\right\}=C_i-{\min }_k\left\{C_k\right\}$$

(6)

This operation removes the same value from all bands: in fact, it amounts to removing the largest possible white component at each pixel while preserving the positivity of all digital counts. In an image taken under natural sunlight, it will help to remove highlights or ‘sun glints’, even though they are not strictly white due to the wavelength dependency of refractive indices and variations in solar radiation [17]. The second step is actually a brightness normalization: the vector $\mathbf{c}$ is divided by the sum of its entries, yielding a vector invariant ${\mathbf{L}}^{(M)}$ with entries:

$$L_i^{(M)}=\frac{c_i}{{\sum }_{j=1}^n{c}_j}=\frac{C_i-{\min }_k\left\{C_k\right\}}{{\sum }_{j=1}^n\left(C_j-{\min }_k\left\{C_k\right\}\right)}$$

(7)

In order to illustrate the effect of these operations, we apply them to the RGB image displayed in Figure 3: it shows a set of plain color wooden toys with different shapes, under natural sunlight. The shapes are fitted so as to cast shadows on each other, each hue is then seen with a great variety of illuminations due to the varied orientations of surface normals as well as casted shadows.

Figure 4 shows both the image resulting from the transformation and the location of transformed pixel values in the RGB space. It can be seen that the spectral invariant ${\mathbf{L}}^{(M)}$ always lies on a triangle-shaped ‘wire-frame’, consisting of three connected rods: the frame is three-dimensional but, just as the color wheel, it basically provides only one degree of freedom. This rather complicated shape results from the choice of substracting the pixel-wise minimum of all bands in the first step: ${\mathbf{L}}^{(M)}$ has always at least one entry equal to zero, but it is not always the same. As a result, it always lies on one of the faces of the unit cube. The normalization step yields another property: the entries sum to 1. If we define the ‘white’ vector $\mathbf{w}={[1,1,\dots ,1]}^T$ having all entries equal to 1, the normalization can be written as a scalar product:

$$\sum _j{L}_j^{(M)}=1⟺{\mathbf{L}}^{(M)}·\mathbf{w}=1=\mathrm{constant}$$

(8)

Hence, ${\mathbf{L}}^{(M)}$ is in a subspace (here, a plane) orthogonal to the ‘white’ direction. If we rotate this ‘white’ axis in order to align it onto the last Cartesian axis ${[0,0,1]}^T$, the triangle lies in the plane ${\mathbb{R}}^2$ and each side can be mapped linearly onto a 120° arc of a circle: the popular RGB color wheel is actually built in the exact same way (Hue–Saturation–Value model).

Obviously, the invariant cannot be computed for gray pixels having the same digital count in all bands, since we would have:

$$L_i^{(M)}=\frac{c_i}{{\sum }_{j=1}^n{c}_j}=\frac{C_i-{\min }_k\left\{C_k\right\}}{{\sum }_{j=1}^n\left(C_j-{\min }_k\left\{C_k\right\}\right)}=\frac{0}{0}$$

(9)

In practice, we remove all ‘nearly gray’ pixels, such that

$$0\le \left({\sum }_{j=1}^n{c}_j\right)<\epsilon$$

(10)

It is worth noting that this will produce gaps in the image.

3.2. Multispectral Hue as a Directional Variable

The definition of the multispectral invariant ${\mathbf{L}}^{(M)}$ given by Montoliu et al. [16] essentially spells out the two basic steps needed to define a multispectral hue. Even though we can interpret these two steps in terms of resaturation and brightness normalization, they are essentially a way of reducing the dimensionality of the problem and minimizing the degrees of freedom for hue analysis. The only drawback of the definition is the topology of the subset in which the invariant lies. Therefore, we propose an alternate definition, which will yield a much smoother topology.

3.2.1. Mathematical Definition

In the resaturation step, we choose to remove the pixel-wise mean of all bands, instead of the minimum. We obtain the vector ${\mathbf{c}}^{\prime }$ with entries:

$$c_i^{\prime }=C_i-\frac{1}{n}\sum _{j=1}^n{C}_j$$

(11)

Then, instead of normalizing with the sum of the resulting entries, we simply normalize this vector with its Euclidean norm, yielding the vector invariant $\mathbf{U}$:

$$\mathbf{U}=\frac{{\mathbf{c}}^{\prime }}{\parallel {\mathbf{c}}^{\prime }\parallel }$$

(12)

The entries of this vector read explicitely:

$$U_i=\frac{c_i^{\prime }}{\sqrt{{\sum }_{j=1}^n{c_j^{\prime }}^2}}=\frac{C_i-\frac{1}{n}{\sum }_{j=1}^n{C}_j}{\sqrt{\left({\sum }_{j=1}^n{C}_j^2\right)-\frac{1}{n}{\left({\sum }_{j=1}^n{C}_j\right)}^2}}=\frac{C_i-{\mathrm{Mean}}_j\left\{C_j\right\}}{\sqrt{n-1}{\mathrm{StDev}}_j\left\{C_j\right\}}$$

(13)

3.2.2. Illustration with the RGB Image

Just as we did for the invariant ${\mathbf{L}}^{(M)}$ of Montoliu et al. [16], in Figure 5 we show the resulting transformed image and the location of transformed pixel values in the invariant space. It is worth noting that the transformation does not preserve the positivity of entries $U_i$. However, since $\mathbf{U}$ has unit norm, it lies in the $[-1,1]\times [-1,1]\times \dots \times [-1,1]$ cube and $\frac{U_i+1}{2}$ will remain inside the positive-valued unit cube; we use these values to reconstruct the R, G, and B components.

This time, it is the resaturation condition that must be interpreted as a scalar product. Centering the values implies:

$$\sum _j{U}_j=0⟺\mathbf{U}·\mathbf{w}=0$$

(14)

Again, the invariant lies in a subspace (a plane in ${\mathbb{R}}^3$) orthogonal to the ‘white’ vector $\mathbf{w}$. Moreover, it has unit norm by virtue of the normalization step: as a result, it directy lies on a great circle (‘tilted’ equator). Rotating the white axis onto the last Cartesian axis will again allow us to discard this last coordinate: we are left with a circle that can be described with only two Cartesian coordinates and only one angular coordinate.

This property will extend in any dimension: the multispectral hue $\mathbf{U}$ will always belong to a subset homeomorphic to ${\mathbb{S}}^{n-2}$ the unit $(n-2)$-sphere embedded in ${\mathbb{R}}^{n-1}$, provided that we perform a final rotation that aligns the ‘white’ direction $\mathbf{w}$ onto the ‘North Pole’ $\mathbf{n}={[0,0,\dots ,0,1]}^T$ of ${\mathbb{R}}^n$. The general expression of the corresponding rotation matrix is given in Appendix A.

3.2.3. Summary of the Section

In this section, we defined a multispectral hue, which is an extension of the classical RGB hue to a higher number of bands: it allows us to resaturate and normalize the brightness of each pixel, whatever the number of bands. [7]. This multispectral hue denoted $\mathbf{U}$ is a direction (unit vector) in a $(n-1)$ dimensional Euclidean space, and it has $(n-2)$ angular degrees of freedom:

• With $n=3$ bands, $\mathbf{U}$ is on ${\mathbb{S}}^1$ the unit circle in the Euclidean plane ${\mathbb{R}}^2$: it is closely related to the ‘color wheel’;
• With $n=4$ bands, $\mathbf{U}$ is on ${\mathbb{S}}^2$ the ‘usual’ unit sphere in the Euclidean space ${\mathbb{R}}^3$. It has two degrees of freedom, we might call them ‘color latitude’ and ‘color longitude’ as will be seen later;
• With $n>4$ bands, the hypersphere ${\mathbb{S}}^{n-2}$ cannot be pictured out. However, it is still a smooth Riemannian manifold, and all familiar properties on ${\mathbb{S}}^1$ and ${\mathbb{S}}^2$ (geodesic distance, Euler angles, rotation group, etc.) can be extended to any dimension.

We can then use directional statistics to analyze hue distributions in any dimension.

3.3. Mixture Models for the Analysis of Multispectral Hue Distribution

The distribution of multispectral hues in a given image can be rather complicated. For this reason, it may be convenient to model the probability density function of this directional variable as a mixture of simpler distributions. On the color wheel ${\mathbb{S}}^1$ for the RGB image displayed in Figure 3, we can use the Von Mises distribution as a basic component; its density reads:

$$f(\phi |\mu ,\kappa )=\frac{e^{\kappa cos(\phi -\mu )}}{2\pi I_0(\kappa )}$$

(15)

where $\phi$ is the hue angle, $\mu$ the mean angle, $\kappa$ a concentration parameter and $I_0$ the modified Bessel function of order 0 for normalization. This density can also be written with vectors; noting

$$\mathbf{U}=\left[\begin{array}{c}cos\phi \\ sin\phi \end{array}\right]{\mathbf{v}}_{\mathbf{m}}=\left[\begin{array}{c}cos\mu \\ sin\mu \end{array}\right]$$

(16)

then

$$f(\mathbf{U}|{\mathbf{v}}_{\mathbf{m}},\kappa )=\frac{e^{\kappa {\mathbf{v}}_{\mathbf{m}}·\mathbf{U}}}{2\pi I_0(\kappa )}$$

(17)

where ${\mathbf{v}}_{\mathbf{m}}$ is now the mean (vector) direction.

3.3.1. Mixture Density

Having chosen the elementary component, the mixture density is simply a weighted sum of m components:

$$f(\phi )=\sum _{j=1}^m{\pi }_j{f}_j(\phi |{\mu }_j,{\kappa }_j)$$

(18)

The resulting density still has to sum to unity on the circle, so the weights must sum to 1:

$${\int }_{\mathbb{S}}f(\mathbf{U})d\mathbb{S}={\int }_0^{2\pi }f(\phi )d\phi =1⟹\sum _{j=1}^m{\pi }_j=1$$

(19)

Figure 6 (left) shows the modeled density of hue in our RGB image: we use a five-component model, one for each wooden shape.

3.3.2. Parameter Estimation Using Expectation–Maximization Algorithm

It is necessary to present briefly how the parameters ${\{{\mu }_j,{\kappa }_j,{\pi }_j\}}_{j=1\dots m}$ of the mixture density are estimated, because we will use a variant of this procedure for depth estimation. The procedure used is the iterative expectation–maximization (EM) algorithm [18]:

• E-step (Expectation): at each iteration $(r)$, we start with current parameter estimates ${\left\{{\mu }_j^{(r)};{\kappa }_j^{(r)};{\pi }_j^{(r)}\right\}}_{j=1\dots m}$ and we define updated membership probabilities using Bayes formula:

  $${\pi }_{ij}^{(r+1)}=\frac{{\pi }_j^{(r)}{f}_j\left({\phi }_i|{\mu }_j^{(r)},{\kappa }_j^{(r)}\right)}{{\sum }_{j^{\prime }=1}^m{\pi }_{j^{\prime }}^{(r)}{f}_{j^{\prime }}\left({\phi }_i|{\mu }_{j^{\prime }}^{(r)},{\kappa }_{j^{\prime }}^{(r)}\right)}$$

  (20)

  where ${\pi }_{ij}^{(r+1)}$ is the estimated probability that pixel i with hue angle ${\phi }_i$ comes from component j, with ${\sum }_{j=1}^m{\pi }_{ij}^{(r+1)}=1$ for each pixel;
• M-step (Maximization): given these membership probabilities, the parameters of each component are re-estimated using a slightly modified maximum likelihood method based on the maximization of:

  $${\ell }_j\left({\mu }_j,{\kappa }_j\right)=\sum _{i=1}^{N_{\mathrm{pix}}}{\pi }_{ij}^{(r+1)}log{f}_j\left({\phi }_i|{\mu }_j,{\kappa }_j\right)$$

  (21)

  We see that the weights in this log-likelihood are the membership probabilities, instead of just 1. ${\ell }_j\left({\mu }_j,{\kappa }_j\right)$ can be maximized for each component j independently, leading to updated parameters $\left({\mu }_j^{(r+1)},{\kappa }_j^{(r+1)}\right)$ and

  $${\pi }_j^{(r+1)}=\sum _{i=1}^{N_{\mathrm{pix}}}{\pi }_{ij}^{(r+1)}$$

  (22)

4. Using Multispectral Hue as a Predictor for Depth

In this section, we test our newly-defined multispectral hue as a predictor for depth on a 35-km reach of the Garonne River between Toulouse and Verdun-sur-Garonne (see Figure 7 left). This reach has a rather natural morphology with few embankments past Toulouse, and many wetlands (it is dubbed the “Garonne débordante”, which means “overflowing Garonne” in French). The average channel width is about 150 m and the maximum depth about 5-6 m at low flow. The channel is mainly bedrock, with typical sculpted molassic forms [19] and several major knickpoints with very shallow flow (depth < 1 m).

This site has been widely used in the French remote sensing community to test discharge estimation algorithms using water surface elevation measurements such as the ones that will be produced by the joint NASA and CNES ‘Surface Water and Ocean Topography’ (SWOT) mission [20,21,22]. As our study contributes to these developments, we needed to produce a detailed bathymetric map on this specific reach, and we present the results in this section. As it will be shown further on, the dataset used has some flaws that prevent us from drawing very general conclusions regarding the advantages of the method. A full qualification of the method is needed and will be the object of a future publication using other publicly available datasets [23].

4.1. Dataset

4.1.1. High-Resolution Airborne Imagery

Our study is based on the BD ORTHO^® database provided by the IGN (French National Geographic Institute). This nationwide dataset consists of RGB and Infrared Color (IRC) ortho-rectified and radiometrically corrected airborne images [24,25] with a resolution of 0.50 m. For our area of interest, we use images acquired during the period 17–30 July, 2016. This dataset provides us with four coregistered spectral bands, namely near-infrared (NIR), red, green, and blue.

4.1.2. In Situ Depth Measurements

As mentioned in Section 2, the relationship between hue and depth can be determined either physically (using radiative modeling) or empirically, using calibration against in situ data. In this study, we use the latter approach: the relationship is calibrated against in situ bathymetric data collected at 120 cross-sections along the reach (about 1 section every 400 m). These data were collected during surveys in the period 1987–2002: hence, this end-of-20th-century bathymetry may differ substantially from present bathymetry locally, even for a bedrock channel (this point will be discussed at the end of the paper).

4.1.3. LiDAR Data in the Floodplain

In order to better constrain the estimation of river depth, we use LiDAR elevation data available in the floodplain. These data come from the BDORTHO^® database also provided by IGN. Even though raster tiles are without gaps, elevation data in the river channel must be ignored: the value at a pixel located in the channel is simply the result of an interpolation between nearby bank pixels. In order to filter out interpolated pixels, the BDORTHO^® elevation rasters come with useful ancillary DST (dist) rasters, which indicate the distance of each pixel to the nearest valid LiDAR point (see Figure 8 right).

4.2. Hue Sphere ${\mathbb{S}}^2$ for $n=4$ Bands

According to the approach presented in Section 3, with $n=4$ the invariant $\mathbf{U}$ does not lie on the circle (1-sphere ${\mathbb{S}}^1$) but on the $(n-2)$-sphere, i.e., the classical sphere ${\mathbb{S}}^2$. In Appendix B, we give the computational detail of the mapping of a 4-band pixel ($\mathbf{C}\in {\mathbb{R}}^4$) onto the unit sphere $\mathbf{U}$ on ${\mathbb{S}}^2$. Here is the final mapping of the four ‘monochromatic’ points, in Cartesian coordinates of ${\mathbb{R}}^3$:

$$\begin{array}{ccc}\hfill {\mathbf{C}}_{\mathrm{NIR}}& =& {[1,0,0,0]}^T⟶{\mathbf{U}}_{\mathrm{NIR}}=\frac{1}{3\sqrt{3}}{[5,-1,-1]}^T\hfill \\ \hfill {\mathbf{C}}_{\mathrm{Red}}& =& {[0,1,0,0]}^T⟶{\mathbf{U}}_{\mathrm{Red}}=\frac{1}{3\sqrt{3}}{[-1,5,-1]}^T\hfill \\ \hfill {\mathbf{C}}_{\mathrm{Green}}& =& {[0,0,1,0]}^T⟶{\mathbf{U}}_{\mathrm{Green}}=\frac{1}{3\sqrt{3}}{[-1,-1,5]}^T\hfill \\ \hfill {\mathbf{C}}_{\mathrm{Blue}}& =& {[0,0,0,1]}^T⟶{\mathbf{U}}_{\mathrm{Blue}}=\frac{1}{3\sqrt{3}}{[-3,-3,-3]}^T\hfill \end{array}$$

In order to illustrate how the hue (hyper)sphere works, we show in Figure 9 where two patches will map, along with the four ‘monochromatic’ points mentioned above. The first patch (A) is composed of water pixels: as absorption is much stronger in IR than in any of the visible bands, it is no surprise that these pixels will map close to the ‘anti-NIR’ hue, corresponding to the multispectral hue of a pixel with digital counts ${[0,1,1,1]}^T$. Conversely, a patch of vegetation (B) will map close to the pure NIR hue, as absorption is not the same in all visible bands but still much stronger than in the NIR.

4.3. Masks

In order to reduce the extent of the area investigated and to narrow the spread of hue distributions, several boolean masks are applied on the raw images in order to retain only river pixels, as illustrated in Figure 10:

• Riparian vegetation as well as floating algae are masked using the classical Normalized Difference Vegetation Index (NDVI) [26], which indicates high chlorophyll content. In practice, the rather low threshold selected (NDVI > −0.3) excludes more than just vegetation pixels: all surfaces with any non-negligible reflectance in the IR will be masked;
• Very dark areas are thoroughly removed: they mainly consist of water pixels in the shadow of riparian vegetation; as these pixels are difficult to characterize spectrally (light has traveled through both leaves and water), the criterion used is simply the mean of the four bands;
• A filter is designed to remove whitewater/white wakes in riffle areas: as these pixels appear in light gray in the visible bands (almost hueless in RGB) but with substantial absorption in the IR (typical digital counts are [0.2 0.6 0.6 0.6]), they are removed on the basis of their low saturation and high brightness in the RGB space along with a threshold in the NIR ($C_{\mathrm{NIR}}<0.3$). Note that they have nothing to do with specular reflections: the surface of water appears white there for any viewing angle, because of the scattering from the bubbles/air pockets, not from surface reflection;
• A mask for man-made features crossing the river, such as bridges or power lines as well as their casted shadows, is built manually (some of these pixels are already masked by the previous filters, but not all).

Once all these masks are combined, a final topological erosion of the resulting mask is performed.

4.4. Hue–Depth Visual Correlation on ${\mathbb{S}}^2$

We now focus on river pixels with a valid multispectral hue. We can first look at the projection onto the sphere of all pixels located at depth measurement points. What we see first in Figure 11 is that they are located close to the anti-NIR hue, as already mentioned above. However, if we plot them with measured depth as color, we see that ‘deep’ pixels tend to concentrate close to a single pole of the sphere. This result is perfectly consistent with Philpot’s Equation (2): as depth increases, the weight of bed reflectance $R_b$ in the mixture reflectance vanishes to zero and the digital counts of river pixels, which are basically reflectances, should tend to the limit vector

$${\mathbf{C}}_{\mathbf{c},\infty }\sim \left[\begin{array}{c}{R}_{c,\infty }({\lambda }_{\mathrm{NIR}})\hfill \\ R_{c,\infty }({\lambda }_{\mathrm{Red}})\hfill \\ R_{c,\infty }({\lambda }_{\mathrm{Green}})\hfill \\ R_{c,\infty }({\lambda }_{\mathrm{Blue}})\hfill \end{array}\right]$$

The hue associated with ${\mathbf{C}}_{\mathbf{c},\infty }$ is thus some point ${\mathbf{U}}_{\mathbf{c},\infty }$ on ${\mathbb{S}}^2$, that we can roughly locate visually. Here, ‘infinitely’ deep means about 5 to 6 m deep.

4.5. Mixture Models on ${\mathbb{S}}^2$ and Higher Dimension, and Depth Predictor

Our last step is to find a quantitative predictor for depth. The two principles are the following:

• We assume that the probability density function of multispectral hue in river pixels can be modeled by a two-component mixture density: the first component will represent the statistical distribution of substrate hue, while the second component will represent the statistical distribution of ‘deep water’ hue;
• We will estimate the parameters of the two components so that the membership probability to ‘deep’ component correlates with depth: the probability ${\pi }_{i,\mathrm{deep}}$ that pixel i with hue ${\mathbf{U}}_{\mathbf{i}}$ comes from the ‘deep’ component will be the predictor for depth $h_i$ at pixel i.

The elementary component used on the unit circle ${\mathbb{S}}^1$ in Section 3 is the Von Mises distribution (‘circular normal distribution’). On the unit sphere ${\mathbb{S}}^2$ and higher dimension, we need an extension of this distribution, similar to a multivariate normal distribution wrapped on the (hyper)sphere. Such an extension is the Fisher–Bingham–Kent (FBK) distribution [27]; for the sake of clarity, basic details about this family of distributions are provided in Appendix C. For $\mathbf{U}\in {\mathbb{S}}^{n-2}$, the density reads:

$$f\left(\mathbf{U}|\left\{{\mathbf{v}}_{\mathbf{j}}\right\},\kappa ,\mathit{\beta }\right)=\frac{1}{\mathcal{C}(\kappa ,\mathit{\beta })}exp\left\{\kappa {\mathbf{v}}_{\mathbf{1}}·\mathbf{U}+{\sum }_{j=2}^{n-1}{\beta }_j{({\mathbf{v}}_{\mathbf{j}}·\mathbf{U})}^2\right\}$$

(23)

The $(n-1)$ vectors $\left\{{\mathbf{v}}_{\mathbf{j}}\right\}$ for $1\le j\le (n-1)$ form an orthonormal basis in which the first vector ${\mathbf{v}}_{\mathbf{1}}$ is the mean direction of the distribution, as f is maximum when ${\mathbf{v}}_{\mathbf{1}}·\mathbf{U}=1$. The normalization constant $\mathcal{C}(\kappa ,\mathit{\beta })$ does not have a closed form, but it can be evaluated numerically [28,29]. The number of degrees of freedom increases almost quadratically with the number of bands ($\sim \frac{1}{2}{n}^2$).

Figure 12 illustrates the decomposition of a given density on ${\mathbb{S}}^2$ into two FBK components. In our case, the two components will be the distributions $f_{\mathrm{bed}}$ of bed substrate hue and $f_{\mathrm{deep}}$ of deep water hue, respectively; for the sake of concision, we denote ${\mathbf{\Theta }}_{\mathbf{bed}}$ and ${\mathbf{\Theta }}_{\mathbf{deep}}$ the parameter set describing each component, so that the mixture density reads:

$$f(\mathbf{U})={\pi }_{\mathrm{deep}}{f}_{\mathrm{deep}}\left(\mathbf{U}|{\mathbf{\Theta }}_{\mathbf{deep}}\right)+(1-{\pi }_{\mathrm{deep}})f_{\mathrm{bed}}\left(\mathbf{U}|{\mathbf{\Theta }}_{\mathbf{bed}}\right)$$

4.6. Modified EM Algorithm for Depth Estimation

The goal of the calibration process is, in the end, to be able to correlate depth $h_i$ at an ‘unmeasured’ pixel i having multispectral hue ${\mathbf{U}}_{\mathbf{i}}$, with the posterior membership probability to the deep component:

$${\pi }_{i,\mathrm{deep}}=\frac{{\pi }_{\mathrm{deep}}{f}_{\mathrm{deep}}\left({\mathbf{U}}_{\mathbf{i}}|{\mathbf{\Theta }}_{\mathbf{deep}}\right)}{{\pi }_{\mathrm{deep}}{f}_{\mathrm{deep}}\left({\mathbf{U}}_{\mathbf{i}}|{\mathbf{\Theta }}_{\mathbf{deep}}\right)\hspace{0.17em}+\hspace{0.17em}\left(1-{\pi }_{\mathrm{deep}}\right)f_{\mathrm{bed}}\left({\mathbf{U}}_{\mathbf{i}}|{\mathbf{\Theta }}_{\mathbf{bed}}\right)}$$

Recall that the ‘posterior’ term means ‘after we get the additional information of pixel hue’. In contrast, ${\pi }_{\mathrm{deep}}$ is the ‘prior’, general probability that a randomly choosen pixel comes from the deep component (it is a constant, corresponding to the weight of the deep component in the mixture distribution): once we obtain the multispectral hue ${\mathbf{U}}_{\mathbf{i}}$ for a particular pixel i, we can update the prior membership probability using Bayes formula to obtain the individual posterior membership probability ${\pi }_{i,\mathrm{deep}}$.

In order for this posterior probability to be correlated with depth, we use a slightly modified version of the expectation–maximization algorithm presented in Section 3, where we first estimated membership probabilities for each pixel (expectation E-step) and then updated the parameters of each component independently using the updated membership probabilities (maximization M-step). Here, we will add an intermediate regression (R-step) between these two steps: instead of feeding the maximization algorithm directly with the membership probabilities ${\pi }_{i,\mathrm{deep}}^{(r+1)}$ for each pixel, we will regress the ${\pi }_{i,\mathrm{deep}}^{(r+1)}$ against measured depth $h_i$ and feed the maximization algorithm with regressed membership probabilities ${\widehat{\pi }}_{i,\mathrm{deep}}^{(r+1)}(h_i)$. Here is the workflow performed at each iteration:

• E-step (Expectation): given current parameter estimates $\left\{{\pi }_{\mathrm{deep}}^{(r)};{\mathbf{\Theta }}_{\mathbf{deep}}^{(r)};{\mathbf{\Theta }}_{\mathbf{bed}}^{(r)}\right\}$ of the mixture at iteration $(r)$, compute the updated probability that pixel i with multispectral hue $U_i$ comes from a deep component:

  $${\pi }_{i,\mathrm{deep}}^{(r+1)}=\frac{{\pi }_{\mathrm{deep}}^{(r)}{f}_{\mathrm{deep}}\left({\mathbf{U}}_{\mathbf{i}}|{\mathbf{\Theta }}_{\mathbf{deep}}^{(r)}\right)}{{\pi }_{\mathrm{deep}}^{(r)}{f}_{\mathrm{deep}}\left({\mathbf{U}}_{\mathbf{i}}|{\mathbf{\Theta }}_{\mathbf{deep}}^{(r)}\right)+\left(1-{\pi }_{\mathrm{deep}}^{(r)}\right)f_{\mathrm{bed}}\left({\mathbf{U}}_{\mathbf{i}}|{\mathbf{\Theta }}_{\mathbf{bed}}^{(r)}\right)}$$
• R-step (Regression): regress deep membership probability as a power law function of measured depth so that

  $${\widehat{\pi }}_{i,\mathrm{deep}}^{(r+1)}=a^{(r+1)}{\left(h_i\right)}^{b^{(r+1)}}$$
• M-step (Maximization): independently maximize the log-likelihood for each component

  $$\left\{\begin{array}{cc}{\displaystyle {\ell }_{\mathrm{deep}}\left({\mathbf{\Theta }}_{\mathbf{deep}}\right)=\sum _{i=1}^{N_{\mathrm{depth}}}{\widehat{\pi }}_{i,\mathrm{deep}}^{(r+1)}log{f}_{\mathrm{deep}}\left({\mathbf{U}}_{\mathbf{i}}|{\mathbf{\Theta }}_{\mathbf{deep}}\right)}\hfill & ⟹\mathrm{updated}{\mathbf{\Theta }}_{\mathbf{deep}}^{(r+1)}\hfill \\ & \\ {\displaystyle {\ell }_{\mathrm{bed}}\left({\mathbf{\Theta }}_{\mathbf{bed}}\right)=\sum _{i=1}^{N_{\mathrm{depth}}}\left(1-{\widehat{\pi }}_{i,\mathrm{deep}}^{(r+1)}\right)log{f}_{\mathrm{bed}}\left({\mathbf{U}}_{\mathbf{i}}|{\mathbf{\Theta }}_{\mathbf{bed}}\right)}\hfill & ⟹\mathrm{updated}{\mathbf{\Theta }}_{\mathbf{bed}}^{(r+1)}\hfill \end{array}\right.$$

  Finally, we update the prior membership probability ${\pi }_{\mathrm{deep}}$ (weight of the deep component in the mixture), which is the mathematical expectation of individual (posterior) membership probabilities:

  $${\pi }_{\mathrm{deep}}^{(r+1)}=\frac{1}{N_{\mathrm{depth}}}\sum _{i=1}^{N_{\mathrm{depth}}}{\widehat{\pi }}_{i,\mathrm{deep}}^{(r+1)}$$

where $N_{\mathrm{depth}}$ is the total number of in situ depth measurement points. Then we iterate back at E-step and so on until convergence is reached.

5. Results

5.1. Error Statistics

The modified expectation–maximization algorithm has several outputs after convergence: the final parameters $\left\{{\pi }_{\mathrm{deep}};{\mathbf{\Theta }}_{\mathbf{deep}};{\mathbf{\Theta }}_{\mathbf{bed}}\right\}$ of the mixture, but also the regression parameters $(a,b)$ such that for each pixel i

$${\widehat{\pi }}_{i,\mathrm{deep}}=a{h}_i^b$$

This output curve is shown in Figure 13, left panel. Taking $h_{\max }=a^{-1/b}$, this relation can be reversed in order to express depth as a function of (known) posterior membership probability:

$${\widehat{h}}_i=h_{\max }{\pi }_{i,\mathrm{deep}}^{1/b}$$

Since ${\pi }_{i,\mathrm{deep}}$ ranges between 0 and 1, estimated depth cannot exceed the ‘saturation’ value $h_{\max }$, which is found to be about 4.9 m. In other words, when hue is close to the maximum of the density of the ‘deep’ component, all we can say is that depth is at least $h_{\max }$, but could be much larger. This is a reasonable behavior: the model will never extrapolate to unbounded values. Figure 13, right panel, shows the resulting scatter plot for estimated vs. measured depth: the RMSE is about 0.59 m and the determination coefficient is $r^2=0.55$.

In Figure 14, we plot the correlation between measured depth and the log of all spectral ratios that can be defined with four bands. Obviously, no single pair of bands exhibits a correlation stronger than the one identified on the hue hypersphere. Most noticeably, none of the spectral ratios is able to identify large depths unambiguously.

A simple benchmark is provided by a Multiple Linear (ML) regression using all independent log-ratios (i.e., ratios between pairs of adjacent bands). It is a specific type of MODPA analysis [15] and the resulting optimal combination of predictors reads:

$$h_{\mathrm{ML}}=2.60log\left(\frac{\mathrm{Red}}{\mathrm{NIR}}\right)+11.48log\left(\frac{\mathrm{Green}}{\mathrm{Red}}\right)-7.31log\left(\frac{\mathrm{Blue}}{\mathrm{Green}}\right)-7.39$$

(24)

The result is shown in Figure 15: the model yields an RMSE of 0.58 m, roughly similar to our procedure. However, the correlation is substantially lower: due to the least-square nature of the model, the variance of estimated depth is lower than the variance of measured depth and the model is not able to yield depths larger than 2.5 m. In contrast, the mixture model on the hypersphere yields depths as large as 4 m.

5.2. Focus on Some Cross-Sections

In Figure 16, we show measured and simulated depth at two representative cross-sections. This illustrates how the method is able to identify localized bed features and high-frequency variations in bathymetry, which are typical of the bedrock channel of the Garonne Toulousaine [19].

6. Discussion

6.1. Sources of Error

6.1.1. Radiometric Inhomogeneity

The BDORTHO^® images used in this study form a nationwide dataset covering all 101 French departments. In practice, these images are radiometrically equalized in order to provide a seamless mosaic at the scale of each department [24,25]. There remain some variations in radiometric aspects inside the mosaic, and especially at the border between two departments. As the investigated reach spans two departments, namely Haute-Garonne (31) and Tarn-et-Garonne (82), the dataset is clearly not radiometrically homogeneous, and this increases the variability in both substrate and deep water hue. This additional variability reduces the precision of the hue—depth relation, but since the method is designed to deal with a dispersion in the distribution of both components through parameter sets ${\mathbf{\Theta }}_{\mathbf{bed}}$ and ${\mathbf{\Theta }}_{\mathbf{deep}}$ of the Fisher–Bingham–Kent distributions, the radiometric inhomogeneity can be statistically treated together with the natural variations. The key parameter acting as a quantification of non-uniformity is the concentration parameter ${\kappa }_{\mathrm{bed}}$ of the bed hue distribution. This is an important advantage of the method over log-ratio methods, which rely on the assumption that reflectance ratios $\frac{R_b(\mathbf{x},{\lambda }_i)}{R_b(\mathbf{x},{\lambda }_j)}$ are rather uniform spatially (see Section 2).

6.1.2. Age and Planar Precision of In Situ Data

As mentioned in Section 4.1.2, the bathymetric data used in this study were collected during surveys in the period 1987–2002. This end-of-twentieth-century bathymetry may differ substantially from present bathymetry locally, even for a bedrock channel. Moreover, a precise positioning of profile with GPS was not always available: errors up to several meters on the transverse position of cross-section points were obvious for a few of them, and likely for many. Profiles were manually repositioned where it was possible, but the remaining uncertainty is an important source of error.

6.1.3. Water Surface Elevation

The airborne images were taken at low flow during the month of July, 2016. However, we do not have a corresponding water surface profile of the reach for the same period; as in situ data consists of absolute bed elevation $z_b$ (in meters a.s.l) and not depth, we needed an estimated map of water surface elevation $z_{ws}$ at low flow, which we derived from a mean of two estimates:

• The first estimate was simply given by the lowest valid elevation at each cross-section in the LiDAR DEM of the floodplain (which also covers the channel through interpolation between the banks). Indeed, the LiDAR survey was conducted during the 2012 low flow period, in hydraulic conditions that we assumed roughly similar to those of the 2016 images;
• The second estimate is the elevation that yields the same width as the one observed in the 2016 images (as defined by the NDWI mask), for each surveyed profile. We then interpolate linearly between the surveyed cross-sections.

We then corrected the values so that our estimated low flow surface elevation profile has strictly monotonically decreasing values. Even if the estimated water surface profile seems reasonable, the uncertainty in water surface elevation, and hence in field estimate of water depth $h=z_{ws}-z_b$, may be several tens of centimeters locally. As large depths drive the parameter estimation, it does not hinder the possibility to identify a robust hue–depth relation, though.

6.2. Perspectives

This work is a proof-of-concept that large-scale identification of river bathymetry is possible using publicly available nationwide datasets. The bathymetric map of the Garonne River produced in this study will be used as an input for reach-scale hydraulic modeling and the analysis of reach-scale hydraulic geometry relationships following [2], with a spatial resolution that could never have been achieved on the sole basis of existing in situ data. However, the shortcomings of the Garonne dataset, listed in Section 6.1, make it necessary to continue the work using higher-quality data.