A Multi-Spatial, Multi-Temporal, Semi-Analytical Model for Bathymetry, Water Turbidity and Bottom Composition using Multispectral Imagery

A reasonable assumption to make for medium to high resolution satellite imagery is constant water turbidity within a spatial region around a given pixel. Klonowski et al [20] used a constant $P$, $G$ and $X$ throughout an entire scene for bottom type classification using hyperspectral imagery.

If we define a spatial region, $r$, around a given pixel, then we (simultaneously) model $N_{r}=(2r+1)^{2}$ pixels. Within this region we constrain $P$, $G$, $X$ and $\Delta$ to be (spatially) static, while $H$, $B_{0}$, $B_{1}$, $\cdots$, $B_{N_{b}-1}$, $q_{0}$, $q_{1}$, $\cdots$, $q_{N_{b}-1}$ may vary.

Thus, we may represent the non-static parameters within each region as

Additionally, in many instances it is reasonable to constrain the variance of H within the spatial region. That is, a constraint of the form $\sigma\left(\textbf{H}\right)<\kappa\,\mkern 1.5mu\overline{\mkern-1.5mu\textbf{H}\mkern-1.5mu}\mkern 1.5mu,$ with $\kappa=0.1$. Where $\sigma\left(\textbf{H}\right)$ and $\mkern 1.5mu\overline{\mkern-1.5mu\textbf{H}\mkern-1.5mu}\mkern 1.5mu$ is the standard deviation and mean of H respectively. (Alternatively, $\kappa$ can be dependent on depth.) The error metric we use in photic for the continuity of H is given by

We can further increase the spectral data available to the model by using a time series of satellite scenes. For LANDSAT-8 and Sentinel-2 scenes, a large number of scenes are available and most are usable for our purposes.

We assume continuity of the seabed over time. Thus, across multiple scenes the depth, $H$, is static up to tidal effects. For $N_{s}$ scenes, it is assumed we know the tidal corrections $H_{\text{tide}_{0}}$, $H_{\text{tide}_{1}}$, $\cdots$, $H_{\text{tide}_{N_{s}-1}}$ to a common datum, such that

We further assume the bottom reflectance is constant over time. This may be a more tenuous assumption than our other assumptions, for example for some coral reefs or seagrass beds.

Thus, the temporal constraints are $H$, $B_{0}$, $B_{1}$, $\cdots$, $B_{N_{b}-1}$, $q_{0}$, $q_{1}$, $\cdots$ and $q_{N_{b}-1}$ are constant, while allowing the water turbidity to vary over time - $P$, $G$, $X$ (and $\Delta$), which we represent as

These constraints complement the spatial constraints, in that previously $H$, $B_{0}$, $B_{1}$, $\cdots$, $B_{N_{b}-1}$, $q_{0}$, $q_{1}$, $\cdots$, $q_{N_{b}-1}$ could vary spatially, while $P$, $G$, $X$ are spatially static.

We will now combine the spatial and temporal extensions into a single error metric for our model. As before, let $N_{s}$ be the number of satellite scenes and $N_{r}$ the size of the spatial region around the spectra we are modelling. Then denote $R_{rs_{{\text{measured}}_{i}}}^{j}(\lambda)$ and $R_{rs_{{\text{modelled}}_{i}}}^{j}(\lambda)$ to be the measured and modelled RRS at spatial location $i$ (for $i=0,\cdots,N_{r}-1$) and scene $j$ (for $j=0,\cdots,N_{s}-1$) for wavelength, $\lambda$, respectively. Then the RMS error across all scenes and the entire region is given by

We have three error metrics - the RMS error (3.4), the SAM error (3.5), and the continuity of H error (3.2). We combine these into a final, weighted, error metric for our model,

where $\omega_{0}\gg\omega_{1}$ and $\omega_{0}+\omega_{1}=1$. In photic, we experimentally found $\omega_{0}=0.85$ and $\omega_{1}=0.15$ balances both terms in the error metric.

In this formulation of the model, we have $N_{r}$ parameters for H; $N_{b}\,N_{r}$ parameters for $\textbf{B}_{0}$, $\textbf{B}_{1}$, $\cdots$, $\textbf{B}_{N_{b}-1}$; $N_{b}\,N_{r}$ parameters for $\textbf{q}_{0}$, $\textbf{q}_{1}$, $\cdots$, $\textbf{q}_{N_{b}-1}$; and $4N_{s}$ parameters for P, G, X, $\mathbf{\Delta}$. In total, we have

parameters to determine.

These spatial and temporal extensions to the model of Lee et al greatly increase data used in the spectral matching. Below we compare the number of measured spectra (assuming we use the LANDSAT 8 satellite with the coastal, blue, green and red spectral bands) across each scene and region to the number of unknown parameters in our model, where the number of bottom types, $N_{b}$ is 2.

Lee et al [15] used $H=10.0$ m as an initial depth estimate and subsequently used an estimate of $H=1/(6P)$ [33]. Starting the model inversion at a reasonable depth estimate was a technique used by Albert et al [22].

If reliable depth profiles are known, then initial, and often accurate, depth estimates can be obtained using an empirical method [4][6][24]. Closely following the model of Lyzenga et al [24] - the empirical depth, $H_{\text{empirical}}$, for a $n$-band spectrum is given by

where $h_{0}=\sum\limits_{i=1}^{n}h_{i}\log(\rho(\lambda_{i-1})/\pi)$ and $\sum\limits_{i=1}^{n}h_{i}\,k(\lambda_{i-1})=1$.

The subsurface RRS of optically deep waters is estimated using the method of Lyzenga et al [24]. We compute both the mean subsurface RRS of optically deep water, $\mkern 1.5mu\overline{\mkern-1.5mur\mkern-1.5mu}\mkern 1.5mu_{rs_{\infty}}(\lambda)$, and the standard deviation of the subsurface RRS of optically deep water, $\sigma\left(r_{rs_{\infty}}(\lambda)\right)$ in each band. It is worth noting that this is a scene-wide estimate.

The attenuation coefficients are estimated using the blue and green spectral bands

then interpolate across spectra and water types from a table of spectral attenuation coefficients for different coastal and oceanic water types to directly estimate the attenuation coefficients and the water type [5][13]. Again, it is worth noting that this is a scene-wide estimate.

Consider $N$ soundings $s_{0},s_{1},\cdots,s_{N-1}$, we begin by generating weights for each soundings $w_{0},w_{1},\cdots,w_{N-1}$ such that

To find the minima of $E_{\text{empirical}}$, we use multiple iterations of the simplex algorithm [2], each with a pseudo-randomly generated initial simplex. This random initialisation with multiple iterations helps prevent the optimisation from settling in a local minimum.

With multiple scenes, the depths from each scene from the empirical algorithm can be synthesised into a single depth estimate by way of a temporal median across all modelled scenes (once we account for tidal variations). Alternatively, a weighted mean can be used. In this case, the weights are inversely proportional to the final $E_{\text{empirical}}$ of each scene.

It would be natural within a computer implementation to model the pixels in a column- or row-major order. However in photic, pixels within the modelling domain are sorted by spectral angle above the deep water mean (via the SAM). The spatial indices of the pixels are stored so we can return the modelling results to their original location. The model starts with the pixels closest to the mean deep water spectrum. We will subsequently see why this ordering is favourable to our modelling.

The initial estimates closely follows Lee et al [15]. Without any knowledge of field samples, the initial parameterisation of P, G, X and $\mathbf{\Delta}$ is given by

where $\mkern 1.5mu\overline{\mkern-1.5muR\mkern-1.5mu}\mkern 1.5mu_{rs}$ is the average over the $N_{r}$ spectra in the region. We call this parameterisation a cold start, as previously modelled pixels with similar spectra have not guided the current starting point.

Alternatively, as the pixels are sorted by their SAM from the deep water mean RRS, if the current region of pixels being modelled is within a model-defined threshold of the previously modelled pixel, then we hot start the model with the previous optimal parameterisation for P, G, X and $\mathbf{\Delta}$. This significantly reduces the number of iterations before convergence.

The initial estimate of $\textbf{B}=0.4{R_{rs_{i}}}(490)$ follows the HOPE model of Lee et al [33]. We arbitrarily set $q_{i}=1.0$, which corresponds to equal proportions of each bottom type. Unlike the initialisation of P, G, X and $\mathbf{\Delta}$, we do not hot start B and q; as this would make assumptions on the continuity of the bottom reflectance.

With our initial estimation of all parameters (both spatially and temporally) in (2.1), we now seek the parameterisation which minimises our error metric, $E_{\text{photic}}$. photic uses the simplex algorithm [2] to perform the optimisation. This algorithm is derivative free and converges quickly if a reasonable starting point (simplex) is chosen.

When the model is cold started and an initial estimate for H is not given, then we perform multiple optimisations with starting points at the following depths 0.1, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.5, 15.0, 17.5, 20.0, 25.0, 30.0 m. This is computationally expensive, but rare as most pixels are hot-started.

photic can also take a range for each parameter. Within the simplex optimisation each parameter is constrained to be within its user-defined range. These ranges could come from local knowledge or from physical measurements from the modelling domain.

Modelling large areas using the optimisation approach requires significant computing resources. One way to decrease this computational burden is the inclusion of a small, dynamic, lookup table (LUT) which is searched prior to running the optimisation approach. If the match between the current pixel and spectra in the LUT (in the sense of the SAM between the two spectra) is below a model-defined threshold then the LUT is used and the optimisation approach is skipped.

Whenever a pixel is modelled with the optimisation approach, if $E_{\text{photic}}$ is below a model-defined threshold then the result of this modelling is stored in the LUT along with a timestamp. When a new entry is stored in the LUT, the oldest entry is discarded (overwritten). In photic the LUT is deliberately small, with only 256 entries. This way searching the LUT is fast, however the LUT should still contain many spectra relevant to the current pixel being modelled as the pixels are sorted prior to modelling.

In photic the threshold on entering a modelled pixel to the LUT is initially $E_{\text{photic}}<\max\left(1.5,1.125\,N_{s}\right)$, but can be adaptively modified. The threshold cannot exceed $2.5+2.5N_{s}$ (where the units of $E_{\text{photic}}$ is relative percentage error). This prevents poorly modelled pixels from further degrading pixels within the modelling domain.

Using a dynamic LUT in conjunction with sorting the spectra results in significant speed improvements. Generally, more than 95% of pixels are modelled using the LUT. The LUT can model in excess of $100\,000\,\text{px}/\text{sec}$, which is at least three orders of magnitude faster than the optimisation approach.

As a detailed example of the model, we give a breakdown of the optimisation process for $N_{s}=2$ and $N_{r}=9$. The spectra are taken from the case study in section 6, for scenes LC81150752018058LGN00 and LC81150752019253LGN00. The location of these spectra are 114.361135E, 21.676824S and from the sounding data the water depth is 14.9 m LAT. The bottom of atmosphere RRS spectra for each scene and region is given in the table below.

All parameters were initialised using the scheme described in sections 4.1, 4.3 and 4.4.

In summary, the model required 1043 iterations to converge to a model error of $E_{\text{photic}}=3.96\%$. The derived parameters were

Thus, the model-derived depth is 16.47 m.

In the first approximation to the inversion, we model all individual scenes, pairs of scenes, 3-tuples of scenes and 4-tuples of scenes. In each of these modelling iterations, all parameters P, G, X, $\textbf{B}_{0}$, $\textbf{B}_{1}$, $\cdots$, $\textbf{B}_{N_{b}-1}$, $\textbf{q}_{0}$, $\textbf{q}_{1}$, $\cdots$, $\textbf{q}_{N_{b}-1}$, H and $\mathbf{\Delta}$ are estimated by the model and stored for subsequent use.

If reliable depth profiles are not known, then we take the median of all iterations of the singletons, pairs, 3- and 4-tuples from the first approximation above. We use a weighted median, where the weight is proportional to the number of scenes in the model. That is, proportional to $N_{s}$.

If reliable depth profiles are known, then the model-derived depths, $H$, can be aligned with these profiles. This technique, for a single scene, was used by Ohlendorf et al [32, Figure 9. Depth validation at Rottnest Island, 2005] as a post processing step to the model inversion.

From the $n$ model-derived depths, $H_{0},H_{1},\cdots,H_{n-1}$, we construct a single depth $H_{\text{aligned}}$, such that

where $c_{i},a_{i},b_{i}>0$. We construct the same weighting scheme as used in the empirical method - consider $N$ soundings $s_{0},s_{1},\cdots,s_{N-1}$, we begin by generating weights for each soundings $w_{0},w_{1},\cdots,w_{N-1}$ such that

where $M=\max(W_{0},W_{1},\cdots,W_{N-1})$. We construct a weighted, root-mean-square, relative error metric, $E_{\text{aligned}}$, such that

where $H_{\text{aligned}_{i}}$ is shorthand for $H_{\text{aligned}}(c_{0},c_{1},\cdots,c_{n-1},a_{0},a_{1},\cdots a_{n-1},b_{0},b_{1},\cdots,b_{n-1})$ at the spatial location of the sounding, $s_{i}$.

To find the minima of $E_{\text{aligned}}$, we use the simplex algorithm [2] with an initial simplex of $c_{i},a_{i},b_{i}=1$ for all $i$. If the model is well calibrated, then we expect this starting point to be close to the minima.

If the modelled attenuation coefficients, $k(\lambda)$ (eqn. 2.2), is relatively constant over time, then we can average $k(\lambda)$ over all model iterations from 5.1. Otherwise all individual model-derived parameters can be averaged.

We perform an additional optimisation iteration to improve the accuracy of the unmixing of multiple bottom types. In this modelling iteration we use the previously estimated P, G, X, $\textbf{B}_{0}$, $\textbf{B}_{1}$, $\cdots$, $\textbf{B}_{N_{b}-1}$, H and $\mathbf{\Delta}$. However, P, G, X and $\mathbf{\Delta}$ are constant. We use the values of $\textbf{B}_{0}$, $\textbf{B}_{1}$, $\cdots$, $\textbf{B}_{N_{b}-1}$ as an initialisation for subsequent modelling.

We begin by expressing 2.1 in terms of $\rho(\lambda)$, which we term $\rho_{\text{modelled}}(\lambda)$ with parameters $P$, $G$, $X$, $H$, $\Delta$

where

and the corresponding unmixed $\rho(\lambda)$, which we term $\rho_{\text{unmixed}}(B_{0},B_{1},\cdots,B_{N_{b}-1},q_{0},q_{1},\cdots,q_{N_{b}-1},\lambda)$ is given by

Then we use the averaged (or aligned) depths as H across all subsequent model iterations. The values of P, G, X, $\mathbf{\Delta}$ are averaged across all previous modelling iterations from 5.1 for each scene. If noise is present in these datasets then we apply a median filter, prior to computing $\rho_{\text{modelled}}$.

In this modelling iteration we keep P, G, X, H, $\mathbf{\Delta}$ static and optimise only for $\textbf{B}_{0}$, $\textbf{B}_{1}$, $\cdots$, $\textbf{B}_{N_{b}-1}$, $\textbf{q}_{0}$, $\textbf{q}_{1}$, $\cdots$, $\textbf{q}_{N_{b}-1}$. As previously, we assume the bottom spectra do not change over time, this allows us to simultaneously model multiple scenes. We index the $N_{s}$ scenes, with $j$.

The error metric for the unmixing is given by

where

and

where $\rho_{\text{unmixed}}^{j}$ is shorthand for $\rho_{\text{unmixed}}^{j}(B_{0},B_{1},\cdots,B_{N_{b}-1},,q_{0},q_{1},\cdots,q_{N_{b}-1},\lambda)$ and $\rho_{\text{modelled}}^{j}$ is shorthand for $\rho_{\text{modelled}}^{j}(P,G,X,H,\Delta,\lambda)$.

Obtaining a realistic depth error estimate requires accounting for uncertainties in both the model inputs and the model itself.

We estimate the sensor noise by computing the RRS standard deviation in optically deep water. The model then adds (or subtracts) random noise within the threshold of the estimated sensor noise. The corresponding standard deviation of the resulting depths over many trials is the depth error estimate.

We can further account for errors in the model, bottom reflectance, absorption and backscattering spectra by scaling the depth error estimate, if reasonable upper bounds on the error estimates of these quantities are known.