Sources and Impacts of Bottom Slope Uncertainty on Estimation of Seafloor Backscatter from Swath Sonars

Abstract

Seafloor backscatter data from multibeam echosounders are now widely used in seafloor characterization studies. Accurate and repeatable measurements are essential for advancing the success of these techniques. This paper explores the impact of uncertainty in our knowledge of the local seafloor slope on the overall accuracy of the backscatter measurement. Amongst the various sources of slope uncertainty studied, the impact of bathymetric uncertainty and scale were identified as the major sources of slope uncertainty. The bottom slope affects two important corrections needed for estimating seafloor backscatter: (1) The insonified area and; (2) the seafloor incidence angle. The impacts of these slope-related uncertainty sources were quantified for a shallow water multibeam survey. The results show that the most significant uncertainty in backscatter data arises when seafloor slope is not accounted for or when low-resolution bathymetry is used to estimate seafloor slope. This effect is enhanced in rough seafloors. A standard method of seafloor slope correction is proposed to achieve repeatable and accurate backscatter results. Additionally, a standard data package, including metadata describing the slope corrections applied, needs to accompany backscatter results and should include details of the slope estimation method and resolution of the bathymetry used.

1. Introduction

Over the last few decades, there have been significant advances in the use of multibeam echo sounders (MBES) to acoustically determine the seafloor properties (e.g., depth, sediment type) that are critical for a host of applications, including nautical charting, seabed habitat assessments, and geological interpretations [1,2,3,4]. While the initial premise of MBES studies has been on extracting bathymetric information, more recently, there is an increasing shift in focus on seafloor backscatter, particularly with respect to those studies aimed at seafloor characterization. Seafloor backscatter studies fall under the general category of remote sensing. With the growing use of remote sensing data, the topic of uncertainty has been receiving increasing attention, particularly in the terrestrial geographical sciences [5,6,7,8]. This focus has led to recommendations that the spatial output of remote sensing data, when compiled into a geographical information system (GIS), should be (at least) twofold: (i) A map of the variable of interest and (ii) some assessment of the measurement uncertainty in that map.

Although calls for including uncertainty estimates in remote sensing studies have been numerous, its practical application is challenging. The main reasons for this are: (1) Knowledge of uncertainty in the measurements is often not available; (2) the impact of the choice of spatial scale is inherently linked to the variable being mapped (which is essentially the unknown), hence the ambiguity as what should be the ideal spatial scale; (3) data processing tools are often not transparent, i.e., their algorithms can be proprietary, prohibiting the calculation of uncertainty propagation from input to the final outputs; (4) the lack of ground-truth data to independently verify the remote sensing observations; and (5) remote sensing data often involves multiple dimensions (at least three in most cases: Position, i.e., x, y, and variable to be mapped). The representation of uncertainty for each pixel thus becomes an issue. Progress on all five challenges is required to improve the uncertainty estimation of remote sensing data. This study is part of a larger effort aimed at looking at various contributors to uncertainty in seafloor backscatter [9] and addresses the uncertainty introduced by uncertain measurements of seafloor slope and the implications of using seafloor slope at different spatial scales for seafloor backscatter corrections.

Two fundamental properties of the seafloor derived from MBES data are the depth and seafloor backscatter strength. A detailed uncertainty model for bathymetry was developed by Hare et al. [10]. Later implementations of uncertainty models, e.g., the combined uncertainty and bathymetry estimator or CUBE [11], have proven that a better understanding of depth estimation and uncertainty can not only improve data quality, but can also simplify data processing. These models identify and capture the bathymetric uncertainty from different sub-components of the sounding measurement, including the echo sounder, the motion sensor, tides, refraction through water, etc. As a result, uncertainty computation methods are now widely adopted across bathymetric data processing tools and end users.

Unlike bathymetric data, which benefits from recognized industry standards used for uncertainty calibration, a “general” uncertainty model for seafloor backscatter has not yet been realized; the sources of uncertainty in backscatter data have only been studied on a case-by-case basis. As techniques to validate backscatter data and treat backscatter quantitatively improve, a need to understand the limitations of backscatter data acquisition and processing has emerged. The major causes of the uncertainty in backscatter data are related to the area correction, incidence angle, transmission loss, and MBES calibration parameters, as well as inherent statistical fluctuations in the measured intensity [9,12]. Two of these corrections (area and incidence angle) depend directly on accurate estimates of seafloor slope. This paper is an extension of an earlier study of the primary factors impacting seafloor backscatter uncertainty [9] and specifically addresses the impact of our ability to measure seafloor slope on the calculation of the area and incident angle with respect to determining seafloor backscatter.

Much of the work on seafloor backscatter uncertainty has been focused on the application of image classification algorithms to backscatter mosaics [13]. For the most part, this work has been conducted without particularly considering the accuracy or quality of the input data. While understanding image classification uncertainty is important and certainly deserves attention, an evaluation of the components of the measurement that create the image data, including geometric and radiometric corrections, must also be included [12,14]. Understanding these physical/geometrical aspects of backscatter measurements and their uncertainties can help improve the backscatter classification and processing. As a community, we need to ensure that we are identifying the steps in backscatter processing that contribute most to uncertainty and focus our attention accordingly. In this context, this paper explores the uncertainties associated with area correction and the incidence angle. Unlike terrestrial studies, where photogrammetric techniques can provide very high-resolution determinations of the terrain slope, the slope derived from bathymetry does not lend itself to direct validation, hence it is difficult to validate an uncertainty model of the seafloor slope. At the same time, it is critical to understand the requirements of the seafloor slope for backscatter corrections and how various choices of seafloor slope computation may impact the seafloor backscatter. This paper addresses some of these questions. Specifically, it: (1) Presents an overview of insonified area corrections and computation of the seafloor incidence angle, and in doing so determines that seafloor slope is the most significant source of uncertainty for these corrections; (2) identifies the major sources of uncertainty for seafloor slope computations; and (3) proposes methods to quantify and reduce the uncertainty of the seafloor slope at appropriate spatial scales in actual real-world surveys.

2. Materials and Methods

MBES backscatter is derived from the measurements of seafloor target strength (see, e.g., [15]). Seafloor target strength is related to the incident and scattered pressure fields from a small patch of the seafloor instantaneously insonified by the sonar signal. The ensemble average (< >) of squared scattered pressure, $<{\left|p_s\right|}^2>$, is proportional to the insonified area, A, and the squared incident pressure, ${\left|p_i\right|}^2$, and inversely proportional to the sonar-target squared distance, $r_s^2$, neglecting absorption and refraction effects:

$$\langle {\left|p_s\right|}^2\rangle ={\left|p_i\right|}^{2}A{\sigma }_b\frac{1}{r_s^2}$$

(1)

where the proportionality coefficient, ${\sigma }_b$, is referred to as the backscattering cross section [16]. Its logarithmic equivalent is the “bottom scattering strength” [15]:

$$S_b=10{\log }_{10}{\sigma }_b.$$

(2)

The target strength (TS in dB re 1 m²) of the seafloor active area, A, is then related to the scattering strength by (for notation simplicity, 10log10A is used instead of the correct form, 10log10 (A/A0), where A0 = 1 m² is the reference unit surface):

$$TS=S_b+10{\log }_{10}A.$$

(3)

In practical situations, the measurements are made using a directional transmitter and receiver. The relation for $S_b$ can be then derived in dB units as [9]:

$$S_b=DN-R{S}_o-S{L}_o+2TL-D_{TX}-D_{RX}-10{\log }_{10}A.$$

(4)

where DN is the received voltage converted and stored as a digital number through an analog-to-digital converter (ADC), $R{S}_o$ is the sensitivity of the receiver transforming the incident acoustic pressure into an electrical signal along its maximum response axis, $S{L}_o$ is the transmit sourced level, 2TL is the two-way transmission loss, $D_{TX}$ and $D_{RX}$ are the transmit (Tx) and receive (Rx) directivity function values in the sonar-target propagation direction, and A the insonified area [17]. For a given seafloor type and frequency, this value of $S_b$ is also related to the seafloor incidence angle, θ.

2.1. Area Insonified Correction

The seafloor patch instantaneously insonified by an MBES signal (defined by the seafloor backscatter sample footprint size) is a function of the transmit and receive beamwidths, as well as a projection of the physical length of the transmitted pulse onto the seafloor. In the classical Mill’s Cross configuration for MBES arrays, the extent of the insonified area in the along-track direction is defined by the Tx sector beamwidth [17]. In the across-track direction, the extent of the insonified area is determined by the Rx (receive) beamwidth in the near-nadir region, and by the pulse length at oblique incidence angles (Figure 1). The accurately detailed computation of the instantaneously insonified area is complicated if both the full Tx and Rx beam patterns are considered; however, approximate formulae are commonly used. This simplification is based on treating the insonified area as the product of the across-track and along-track extent of the projection of the backscatter sample footprint on the seafloor. This approximation for the instantaneously insonified area has been shown to create negligible uncertainty except in the near-nadir region and for wide beams [9,14]. With most MBES operating with beamwidths <1 to 3°, this approximation is not considered a significant source of uncertainty. At an oblique incidence (short-pulse regime, see [17]), the insonified area can be approximated as:

$$A\approx \phi R\frac{cT}{2\mathit{sin}\theta \mathit{cos}{\beta }_y}$$

(5)

and around normal incidence (long-pulse regime) as:

$$A\approx \phi \omega R^2\frac{1}{\mathit{cos}\theta \mathit{cos}{\beta }_y}$$

(6)

where R is the range; φ and $\omega$ are the along-track and across-track two-way equivalent apertures, respectively (see Figure 1); T is the effective pulse length; c is the local sound speed; θ ισ the across-track incidence angle; and ${\beta }_y$ is the along-track seafloor slope. The pulse length, T, considered here, is the equivalent length of either the physical pulse duration in case of continuous waves (CW) or the compressed pulse duration after matched filtering in the case of frequency modulated (FM) transmitted signals [17]. It is to be noted that sonar systems currently do not correct for the seafloor slope in real-time. Instead, a flat horizontal seafloor is assumed during data acquisition [18]. Therefore, the accurate estimation of the insonified area falls on the end-users either by applying third-party software during the post-processing phase or by developing customized software themselves.

2.2. Seafloor Incidence Angle Correction

The seafloor incidence angle is defined as the angle between the receiving beam direction and the vector normal to the seafloor insonified patch while the grazing angle is defined as the angle between the beam direction and the vector parallel to the seafloor patch. The normal vector to the surface patch can be defined as [19]:

$$\overrightarrow{n}={[\tan {\beta }_x,\tan {\beta }_y,-1]}^T$$

(7)

where ${\beta }_x$ and ${\beta }_y$ represent seafloor slopes in the across-track (x) and along-track (y) directions, respectively, and T denotes the transpose operation. The seafloor slope at the beam footprint is computed from the bathymetry in the across-track and along-track directions using either the soundings from the individual beams or a gridded data set. Assuming a flat seafloor, the equation for the nominal receiving beam direction is:

$$\overrightarrow{m}=\left[0,-\cos {\theta }_g,\sin {\theta }_g{]}^T=\right[0,-\sin \varphi ,\cos \varphi {]}^T$$

(8)

where $\varphi$ = nominal transmission angle, and ${\theta }_g=\frac{\pi }{2}-\varphi$ is the nominal seafloor grazing angle. The true incidence angle is then the angle between the two vectors, $\overrightarrow{n}$ and $\overrightarrow{m}$, and can then be given by:

$$\cos {\theta }_{inc}=\frac{\overrightarrow{n}·\overrightarrow{m}}{\Vert \overrightarrow{n}\Vert \Vert \overrightarrow{m}\Vert }$$

(9)

where $·$ is the scalar product and ‖ ‖ is the norm of the vector. The nominal transmission angle (with respect to the sonar) and seafloor incidence angle (${\theta }_{inc}$) are nominally complementary angles unless impacted by refraction through the water column and vessel motion. Hare et al. [10] showed that beam pointing errors at the sonar head depend on uncertainty in a ship’s motion and an estimation of the sound speed at the sonar head. Both of these variables have been well modeled with respect to bathymetric data uncertainty [10] and will not be discussed here as their magnitude as well as their effect on the seafloor incidence angle is usually small.

2.3. Estimation of Seafloor Slope and Its Uncertainty

The seafloor slope “between points” or “for a patch of seafloor” is defined as the rate of change of depth over distance. Considering $z_1$ and $z_2$ as vertical coordinates and $x_1$ and $x_2$ as horizontal coordinates, the slope ($g$) can be calculated as:

$$g={\tan }^{-1}\left(\frac{z_1-z_2}{x_1-x_2}\right).$$

(10)

Applying the law of propagation of variance, the variance in slope is given by:

$${\sigma }_g^2={\left(\frac{{\partial }_g}{{\partial }_{x_1}}\right)}^2{\sigma }_{x_1}^2+{\left(\frac{{\partial }_g}{{\partial }_{x_2}}\right)}^2{\sigma }_{x_2}^2+{\left(\frac{{\partial }_g}{{\partial }_{z_1}}\right)}^2{\sigma }_{z_1}^2+{\left(\frac{{\partial }_g}{{\partial }_{z_2}}\right)}^2{\sigma }_{z_2}^2$$

(11)

where ${\sigma }_X$ represents the respective error in a quantity, X. Hence, estimates of vertical sounding uncertainty (${\sigma }_{z_1}, {\sigma }_{z_2})$ and position uncertainty (${\sigma }_{x_1},{\sigma }_{x_2}$) are required to estimate slope uncertainty.

2.3.1. Directional Slope (Along-Track and Across-Track from Reference Grid)

For backscatter corrections, the two directional slopes (i.e., along- and across-track) are required. As a preliminary estimation, the depth values obtained within each beam can be used to calculate across-track slopes. This calculation is straightforward to implement in the framework of MBES sounding data processing. A frequently-used option is to take at least two depth measurements on either side of a depth sounding and then fit a plane through the depth values. Unedited individual soundings may have large variance that will adversely affect the slope estimation. Therefore, it is important to edit the obvious outliers in the bathymetry before attempting the across-track slope calculation.

While using beam depth data is an option for MBES systems that provide co-located bathymetry and backscatter, the common choice for computing slope is to use the gridded bathymetry. To determine the across-track slope, the gradients along two principal axes, x and y, are calculated first. From the values of dz/dx and dz/dy, one can then determine the slope of the plane that approximates the surface at the local depth point (i, j) using the formula [20]:

$$g_{dir}\left(\gamma \right)={\tan }^{-1}\left[\frac{dz}{dx}\sin \gamma +\frac{dz}{dy}\cos \gamma \right]$$

(12)

where $\gamma$ represents the across-track orientation of the MBES swath, determined based on the vessel heading, along which the slope was calculated. Directional slopes can be computed by applying a 3 × 3 kernel indexed to each cell (except those at the grid edges). Figure 2 illustrates the adopted convention for cell indexing.

Various methods are available for the estimation of directional gradients (dz/dx and dz/dy) from gridded data sets. For a comparison of the performance of these methods, see [21]. Two commonly recommended methods are used here: The central difference method [20] and the Horn [22] method. In the case of the central difference method, the values of two neighboring cells are used to compute gradients along the two axes (Figure 2):

$$\begin{array}{l}\frac{dz}{dx}=\left[z_E-z_W\right]/2S\\ \frac{dz}{dy}=\left[z_N-z_S\right]/2S\end{array}$$

(13)

where S defines the grid cell size. In the case of cells at the grid edge, the central cell itself is used to provide the value for the missing data. For the Horn method, values from six neighboring cells are utilized to compute gradients (Figure 2). The weight applied to each cell depends on its position relative to the central cell. Formulae for the computation of directional gradients in this case are:

$$\frac{dz}{dx}=\frac{\left[(z_{NE}+2z_E+z_{SE})-\left(z_{NW}+2z_W+z_{SW}\right)\right]}{8S}$$

(14)

$$\frac{dz}{dy}=\frac{\left[\left(z_{NE}+2z_N+z_{NW}\right)-\left(z_{SE}+2z_S+z_{SW}\right)\right]}{8S}.$$

(15)

2.3.2. Estimation of Slope Uncertainty

Studying uncertainty in seafloor slope estimation is relatively difficult as the true value of the slope is hard to determine and validate. Dolan and Lucieer [21] and Lucieer et al. [23] recommended using a Monte Carlo simulation to analyze the uncertainty in a slope by considering the uncertainty in the bathymetry. Therefore, the uncertainty of bathymetry data estimated using the CUBE algorithm [11] was used to infer uncertainty in an estimated slope. Besides the horizontal and vertical uncertainty of the soundings, the slope estimation may also be affected by the method used to compute the seafloor slope. Several different algorithms are available for slope estimation that differ based on the method used to determine the neighboring grid cell sizes [24], resulting in different estimates. Additionally, the scale of the seafloor slope (i.e., at what spatial resolution the seafloor slope is estimated) has been shown as an additional source of uncertainty [25] if the appropriate spatial scale is not chosen for a given application.

Although the uncertainties in the soundings and resulting depth grids can be estimated based on bathymetry uncertainty models, their propagation to the slope values is not straightforward. Consider a few soundings in Figure 3, where an arbitrarily assumed horizontal and vertical uncertainty has been used to plot the error ellipses. The slope between any two points is then the slope of the line joining the two error ellipses. As there can be infinite cases for these lines, uncertainty propagation using the Taylor series expansion method is not feasible and use of the Monte Carlo method has been recommended [23]. Therefore, such a simulation was used during this study to estimate the effect of bathymetric uncertainty and the slope computation algorithm on slope estimation. For depth grids, the position uncertainty was assumed as 0, as the depth grid node locations are fixed, and only the vertical uncertainty of the grids was considered. For each simulation iteration, a new bathymetric grid was assembled by adding randomly drawn bathymetric uncertainty values to the mean depth values at each grid node.

2.4. Multibeam Sonar Test Dataset

To illustrate the uncertainty introduced by the processing methodology described above, data from an EM 3002 dual-head MBES collected in water depths of 0.5 to 25 m (in vicinity of Portsmouth, NH) were used. The EM 3002 [26] is a shallow-water multibeam system operating at a nominal frequency of 300 kHz. It transmits a pulse of 150 $\mathsf{\mu }$s and nominally forms 160 beams. The different methods of seafloor slope estimation were applied to the test data and their differences were noted.

MBES Bathymetric Data Processing

The bathymetric data were processed using QIMERA [27]. Tide data were applied based on the Hampton Road, NH tide gauge located ~6 km from the survey area. Automated data editing and cleaning were performed using CUBE [11], which provided estimates of horizontal and vertical total propagated uncertainty. The data were successively gridded at a 0.5, 1, 2, 5, 10, and 20 m grid resolution. These six bathymetric grids represent the best estimate of depth at the respective grid scales. CUBE also generated an uncertainty layer for the bathymetric grids, at different resolutions, representing one standard deviation of the depth uncertainty at each grid node. An overview of the bathymetry of the study area is shown in Figure 4. The survey area consisted of rock outcrops in the northern section of the survey, while the southern section was primarily a flat sandy area. For backscatter data processing, the Fledermaus Geocoder Toolbox (FMGT) [28] was used. FMGT was selected for backscatter processing as it can incorporate user defined bathymetric grids for area and incidence angle corrections, thus allowing a study of the impacts of the choice of grid cell size.

3. Results

Uncertainty of the insonified area and incidence angle caused by the seafloor slope uncertainty is scaled with survey depth and sonar characteristics. To illustrate the uncertainty caused by the seafloor slope and to help quantify the impact of seafloor slope on these corrections, a few examples are presented below based on typical depths where the EM 3002 is operated.

3.1. Slope Impact on the Footprint Extent

For a given beam aperture and pulse length, the area correction increases with depth and the resulting area correction can vary in the range of −10 to 15 dB re.1 m² (Figure 5).

The across-track and along-track slopes, if present, can cause significant variations in the above estimation of the insonified area correction and therefore cannot be ignored when computing the insonified area. Examples of the expected effects of the across-track slope (Figure 6) and along-track slope (Figure 7) are provided for a depth of 100 m. Both figures show that an across-track slope has the more significant effect on the insonified area extent in comparison with the along-track seafloor slope. This could be expected from Equations (5) and (6): The effect of the along-track seafloor slope appears as $\cos {\beta }_y$, hence causing modest variations (0.62 dB for an extreme slope value of ${\beta }_y={30}^{°}$).

3.2. Seafloor Slope Impact on Incidence Angle

If not corrected for the seafloor slope, the seafloor incidence angles will be subject to errors corresponding to the slope magnitude. Assuming the seafloor slope is ignored in the estimate of the incidence angle, the along-track slope is shown to have a greater effect on the near-nadir beams than on the outer beams (Figure 8), while as expected, the across-track slope is shown to have a similar effect on all beam angles across the swath (Figure 9).

The above results clearly illustrate that the seafloor slope is an important factor in the estimation of the insonified area as well as of the incidence angle. The next sections will focus on the results obtained with differences in the scale followed by the uncertainty observed due to the use of different slope estimation methods and vertical uncertainty in soundings.

3.3. Scale Dependent Slope Estimation Uncertainty

The Horn method [22], as implemented within ARCGIS [29], was used to compute the seafloor slopes for bathymetric grids at various grid cell size resolutions, on the data presented above in Figure 4. The differences in the features that can be resolved in the slope layer as a function of the resolution are apparent in Figure 10, with a larger grid cell size smoothing out the detailed seafloor topography. The effect is most prominent in the rough area.

Increasing the grid cell size (from 1 m to 20 m) results in a decrease in the spread (range and standard deviation) of the slope values (

Table 1. Comparison of seafloor slope statistics obtained for different grid cell sizes (see Figure 10). MIN refers to the minimum slope, MAX refers to the maximum slope, RANGE show the differences between MIN and MAX while MEAN is the average and STD is the standard deviation in each cell.

	Grid Cell Size	MIN	MAX	RANGE	MEAN	STD
Flat area	20	0.02	1.02	0.99	0.24	0.14
Flat area	15	0.00	1.22	1.22	0.26	0.16
Flat area	10	0.01	1.84	1.84	0.28	0.18
Flat area	5	0.00	3.44	3.44	0.35	0.29
Flat area	1	0.00	16.70	16.70	0.78	0.64
Rough area	20	0.06	6.48	6.42	2.10	1.30
Rough area	15	0.05	7.33	7.28	2.24	1.33
Rough area	10	0.02	17.24	17.22	3.32	2.14
Rough area	5	0.03	29.98	29.95	4.51	3.18
Rough area	1	0.00	45.07	45.07	5.74	4.39

). This is intuitive since enlarging the grid cell size effectively acts as a low-pass spatial filter on the depth; the standard deviation of the slope estimates is hence lowered, but the information content of the slope estimates is also reduced. The uncertainty due to the choice of the cell size therefore cannot be directly estimated from seafloor slope statistics at a single scale. This is an important realization and emphasizes the fact that during the seafloor backscatter processing, the scale of the seafloor slope estimation is a key factor.

Recognizing that the grid cell size has a strong impact on the slope values, the question of an optimal choice of the grid cell size for the computation of the seafloor slope arises. This is a fundamental issue, which essentially depends on the specific application. Instead, this question can also be approached practically from the perspective of the highest resolution possible from a given data set. In regard to backscatter corrections, the seafloor slope used to determine the incidence angle and insonified area should ideally be at a scale comparable to the backscatter sample footprint. Near nadir, the footprint is defined by the combined Tx-Rx beamwidth for both the across- and along-track extent. Away from nadir, the across-track extent is fixed and depends on the pulse length projection in the across-track direction, while the along-track extent is controlled by the received beamwidth and depth. However, in contrast to backscatter samples, the soundings are usually obtained by considering several adjacent time samples, resulting in a bathymetric resolution that is typically less than the resolution of the backscatter. Subsequently, the best bathymetric resolution that can be obtained from the depth data is actually limited by the beam spacing. Moreover, to compute a bathymetric grid, the depth points are averaged together, resulting in low noise, but also lower resolution. Therefore, the bathymetric data derived either from individual soundings or from grids (even constructed at the highest possible resolution) may be inadequate for the slope correction of the backscatter data. Additionally, the choice of larger grid cell sizes may be necessary to simplify processing and to reduce noise in the bathymetry-derived slopes.

3.4. Uncertainty Due to the Slope Estimation Method

An assessment is made here of two different methods for estimating slopes using either the bathymetric depths within the measurement swath for each ping, or the gridded bathymetry. For gridded bathymetry with a 1 m cell size, the slope was evaluated using two different computing schemas, the central difference method and the Horn method. The beam bathymetry will invariably provide a better resolution as compared to the gridded data set. The expected resolution from beam bathymetry depends on the depth and beam spacing defined by the MBES settings. Commonly used beam-spacing modes include equiangular or equidistant placement of soundings on the seafloor. Typical beam-spacings for various depths for near nadir and outer beams are given in

Table 2. Across-track beam spacing and backscatter sample footprint for various depths and transmission angle for MBES with a beamwidth of 1.5° and pulse length of 150$\mathsf{\mu }\mathrm{s}$ assuming a sound speed of 1500 ms⁻¹. The across-track backscatter footprint depends on the depth in near nadir region, but only depends on the pulse length at oblique angles.

Across Track Beam Spacing (m)					Across Track Backscatter Footprint
	Equiangular				Equidistant
Depth (m) ${\theta }_{inc}$	0°	45°	60°	All angles	Near Nadir	Oblique Angles
10	0.26	0.53	1.09	0.26	0.26	0.11
20	0.52	1.06	2.19	0.53	0.52	0.11
50	1.30	2.68	5.48	1.34	1.30	0.11
100	2.61	5.37	10.97	2.68	2.61	0.11
200	5.23	10.75	21.94	5.36	5.23	0.11

for an MBES with a beamwidth of 1.5° and 160 beams using the equiangular and equidistant mode for a swath of ±65°. As can be seen in Table 2, the highest possible resolution of bathymetry obtained from beam bathymetry is insufficient for backscatter corrections.

Using the previously mentioned MBES survey data in water depths of ~10 m, large differences (as much as 15°) are apparent between the slopes estimated from the beam bathymetry vs. the slopes derived from gridded bathymetry (Figure 11) for a single ping chosen in the rough area, whereas the differences between the two grid-based methods are comparatively small (Figure 12). The differences in the slopes from the beam bathymetry and gridded data set were comparable for the flat areas (not shown). This observation indicates that the bathymetry type (individual soundings vs. gridded bathymetry) has a more prominent effect on the slope than the choice of the algorithm used to estimate the gridded slope. Considering that bathymetric detail is lost when using bathymetric grids (especially with larger grid cells), it is recommended that individual soundings be used when computing seafloor slopes for backscatter corrections. It is also to be noted that individual soundings were not edited for outliers for this comparison. This did not result in spikes in the estimated slopes in this case, but this can be a major issue when excessive noise is present in the individual soundings. The large variability in the slopes from individual soundings is not an indication of large uncertainty, but it more closely captures seafloor slope changes. When the bathymetric soundings are excessively noisy, the only choice is to use the highest possible resolution of the bathymetric grid.

The difference between the slope computation algorithms was minimal in flat areas, but reached ~1° for rough areas. Therefore, the choice of algorithm is considered a minor source of slope uncertainty, but cannot be totally ignored.

3.5. Propagation of Depth Uncertainty to Slope Uncertainty

Slope values are directly impacted by the horizontal and vertical uncertainty of depth measurements. The Horn method described in the previous section is used here to illustrate the propagation of depth uncertainties to slope values. As the horizontal positions of the grid nodes are fixed at regular intervals, the horizontal component is not included in the uncertainty propagation. The CUBE algorithm implemented in QIMERA was used to derive the vertical uncertainty for each grid node at a grid cell size of 1 m (Figure 13).

The slope for each grid point was computed over 50 iterations. For each iteration, the vertical uncertainty values were randomly drawn from the range of CUBE-generated uncertainties and a new bathymetric grid was assembled. The resulting grid was then used as input to compute the seafloor slope using the Horn method.

For the example provided in Figure 14, the slope standard deviation is related to the slope value of the seafloor: For steep seafloors, the value of standard deviation is larger (2–3°), and for flat seafloors, the slope standard deviation is comparatively small (<2°). It is to be noted that uncertainty due to the bathymetric uncertainty is related to the seafloor slope as well. The areas with rough seafloor are expected to exhibit higher bathymetric uncertainty due to large slope variations. The Monte Carlo simulation has been shown as a potential approach to compute seafloor slope uncertainty due to bathymetric uncertainty [21]. However, estimates of slope uncertainty due to bathymetry uncertainty cannot be generalized, and for each data set, the simulations (e.g., using the Monte Carlo method) will need to be carried out to compute slope uncertainty. This poses a challenge for end users who may not have resources to conduct such simulations.

3.6. Impact of Unresolved Seafloor Slope on Backscatter Ensemble Average

When considering the impact of low-resolution bathymetry on backscatter uncertainty, the effect of the scale of the bathymetry on ensemble averaging of the backscatter must be considered. To avoid bias in the backscatter strength, the averaging needs to be at a scale where the individual samples can be assumed to be derived from the same seafloor patch. Currently, a widely-accepted standard does not exist to define the bin size for ensemble averaging and a large range of bin sizes are used: Ranging from angular bins of 0.5° to 1°, with cell sizes much smaller than the corresponding bathymetry for the creation of high-resolution mosaics. The slope-related insonified area and incidence angle uncertainties for individual backscatter samples are therefore propagated to the backscatter strength estimate, with their magnitude depending on the variance of the unresolved seafloor slope. The variance of unresolved slopes, however, is a function of local topography and cannot be determined using the MBES bathymetry. Despite the inability to demonstrate this empirically, it is expected that the bathymetry (at a particular scale) enables computation of the mean slope reliably and the variations around the mean slope will likely follow a normal distribution. Thus, the spatial scale of the slope should practically provide the lower limit of the scale, below which ensemble averages should not be carried out. For instance, if one samples the slope values at the scale of the bathymetry grid, one cannot build $S_b \left(\theta \right)$ by using backscatter elementary values at a finer scale, where the seafloor slope is not resolved.

Although not known, but assuming an arbitrary value of variance in unresolved seafloor slopes, a first-order estimate of the resulting uncertainty in backscatter can be made. To demonstrate this numerically, the values of modeled backscatter have been calculated using the Generic Seafloor Acoustic Backscatter (GSAB) model [30] and used to estimate the effect of angular averaging. In its simplest form (a Gaussian law for specular regime and Lambert’s law at oblique incidences), the impact of averaging over $\pm \theta$ is plotted in Figure 15. The expression for ${\sigma }_b$ for the GSAB model is given by:

$${\sigma }_b=A\exp \left(-\frac{{\theta }^2}{2B^2}\right)+C{\cos }^D\theta$$

(16)

where A is the specular maximum amplitude, B is the facet slope standard deviation, C quantifies the average backscatter level at an oblique incidence, and D is the backscatter angular decrement [30]. The two cases illustrated in Figure 15 are typical angular backscatter curves for (a) a soft-sediment with high narrow specular and oblique decrease in ${\cos }^2\theta$, and (b) a coarse sediment with a low and wide specular and oblique decrease in $\cos \theta$. For most seafloors, the oblique-regime average angle dependence lies between the cosθ and cos²θ curves shown here. The model input parameters (A, B, C, D) are, respectively, (0.1; 2°; 0.001; 2) and (0.03; 7°; 0.01; 1). The values of backscatter obtained from the GSAB model were first computed for each incidence angle, then binned (at various angular bins of 1° to 10° corresponding to the assumed variance of the unresolved slope) and reported at the central angle of the bin. As expected, the resulting uncertainty is maximum for the specular regime, where its magnitude depends on the gradient of the specular lobe. On the other hand, the sensitivity to unresolved seafloor slopes becomes negligible in the “plateau” regime (>15° in this example) and at larger incidence angles. This implies that averaging backscatter at a scale larger than the backscatter sample footprint (e.g., at the scale of gridded bathymetry cells) does not cause a detectable bias even if the small-scale relief causes significant slope variations (as high as 10°) inside the grid cells around the average slopes determined from bathymetry.

On the other hand, the exact impact of unresolved slopes on the insonified area correction cannot be predictively modeled using GSAB. However, if the highest possible resolution of bathymetry is used (i.e., the individual soundings), the variations in the slope are expected to result in negligible bias in the computed ensemble average over random variations in the insonified area. Hence, averaging backscatter at the scale of the gridded bathymetry is feasible except for the specular region (Figure 15a), as it provides bias-free smoothed backscatter at a resolution consistent with bathymetry.

3.7. Practical Impact of Slope Scale on Incidence Angle and Processed Backscatter Results

The MBES survey data processed using different spatial scales for the seafloor slope correction showed that the incidence angle could vary up to tens of degrees depending on the spatial scale of the seafloor used to estimate the incidence angle. To look at the impact of the slope scale on the incidence angle, QPS FMGT was used as it enables the use of a reference grid to calculate incidence angles. For comparison, the results were exported as text files using the built-in export functionality of QPS FMGT (‘ASCII ARA beam detail’). The exported data included the true incidence angle as well as corrected backscatter data after applying corrections for the seafloor slope (based on whether a reference grid was used) and insonified area. For both incidence angles and backscatter levels, the largest differences were observed between no-slope corrections and slopes computed at a 1 m grid resolution (Figure 16). The backscatter values showed variations of up to 3 dB around the specular area, decreasing at oblique angles (Figure 16). As QPS FMGT currently does not provide details of area corrections separately, these differences show the cumulative effect of the backscatter corrections, including the insonified area. However, as the only difference in processing was the change in the reference grid resolution, the differences in the backscatter results imply that these are related to insonified area corrections.

4. Discussion

The focus of this research was on understanding the sources of slope uncertainty while correcting seafloor backscatter data for an insonified area and incidence angle. This cannot be separated from the issues of the analysis scale and choice of methodology for creating underlying bathymetry that is used to correct the backscatter data for the slope and incident angle. As multibeam survey acquisition becomes prevalent, the use of multiple backscatter data sets from different sources is also becoming more widespread. The compilation of backscatter data from various sources, collected using multiple MBES systems, and processed with different tools, highlights a need to understand and quantify uncertainty in the backscatter data. The results presented here have demonstrated that the slopes estimated for correcting seafloor backscatter vary depending on the computational algorithm and uncertainty of the bathymetry. The greatest uncertainty is introduced if the seafloor slope is not resolved at an appropriate scale. Unfortunately, documentation accompanying backscatter results often lacks a description of how the seafloor slope was corrected for, making it difficult for end-users to determine uncertainty due to slope corrections. Information about the computation algorithm and bathymetric data used, and the scale of analysis should accompany seafloor backscatter data so that users can interpret the backscatter data and their intrinsic uncertainty.

4.1. Summary of Uncertainty Components of Seafloor Backscatter Measurements Related to Seafloor Slope

The impact scale adapted from [9] was used to classify the magnitude of the uncertainty related to area correction (in dB) and incidence angle (in °) (

Table 3. Scale adapted from [9] to classify the magnitude of uncertainty in area correction and incidence angle.

	Negligible (N)	Small (S)	Moderate (M)	High (H)	Prohibitive (P)
Area correction (dB)	0.01 to 0.1	0.1 to 1	1 to 3	3 to 6	Beyond 6
Incidence angle (°)	0.01 to 0.1	0.1 to 1	1 to 3	3 to 6	Beyond 6

). The impacts of the main uncertainty sources in seafloor slope computation are summarized in

Table 4. Major sources of uncertainty for seafloor slope required for area insonified and seafloor incidence angle. See the code (N-S-M-H-P) definition in Table 3.

Seafloor Slope Uncertainty Source	Magnitude of Uncertainty in Area Insonified and Incidence Angle	Possible Quality Improvement
Flat seafloor assumption (seafloor slope completely ignored)	N to P depending on topography	Use bathymetry in slope compensation
Inappropriate scale of seafloor slope computation (beam bathymetry vs. grid at lower resolution)	M to H depending on large scale topography	Use of highest available resolution bathymetry
Unresolved seafloor slope	N to M depending on small scale topography	Average backscatter values inside angular bins
Bathymetry uncertainty	S to M depending on bathymetric uncertainty and magnitude of seafloor slope	-
Seafloor slope algorithm based on bathymetry grid	N to S	-

.

4.1.1. Approaches to Using Bathymetry for Slope Estimation

As is evident from the discussion above, the largest uncertainty levels are obtained when the seafloor slope is ignored while processing backscatter data. Several works have highlighted this issue [19]. This uncertainty directly depends on the slope magnitude (i.e., on the local topography); therefore, its prediction cannot be accomplished a priori. However, as the collection of concurrent bathymetry is conducted during MBES backscatter surveys, this uncertainty source can be reduced if the bathymetry data are used to compute the seafloor slopes. At this point, the issue of the bathymetry spatial scale must be carefully considered. Terrestrial studies of slope accuracy support the conclusion that spatial scale is an important factor for accurately estimating the terrain slope [25,31], which has also been shown to be critical for seafloor slope [21]. Several options, with varying resolutions, are available for end-users to compute the slope from the available bathymetric data. An important question is whether to use across-track beam bathymetry (pros: Optimal resolution; cons: Noise unfiltered; does not give the across/along-track slope at consistent spacing) or gridded bathymetry (pros: Gives along and across-track slopes with a consistent sampling spacing, and well filtered; cons: The gridding is often too coarse to provide detailed slopes). The results from Figure 7 indicate that along-track slope effects are relatively small (0.1–1 dB) for an area that is insonified and affects only the near-nadir beam adversely. In contrast, the results from Figure 11 indicate that estimated slope values vary significantly depending on whether beam bathymetry or a bathymetric grid is chosen for the slope estimation, especially for seafloors with large topographic variations. The uncertainty in the area insonified and incidence angle can be prohibitively high (> 6 dB; > 6°) as indicated in Figure 16 for complex terrains with rough topography. Therefore, for the across-track slopes, beam bathymetry is the preferred choice. For the along-track slope, the gridded bathymetry at the highest resolution can offer a reasonable solution. Hybrid processing of bathymetry (using beam bathymetry for the across-track slope and using gridded bathymetry for the along-track) could also provide an optimal solution.

4.1.2. Impact of Spatial Scale

For the insonified area computation required for backscatter processing, the spatial scale at which the slope needs to be resolved should be ideally the backscatter sample footprint extent. The interplay between calculable uncertainties vs. lack of information caused by using a larger grid cell size will remain an unanswered question as generalized terrain characteristics cannot be reliably assumed for a surveyed area. However, it is reasonable and useful to estimate upper bounds for this cause of uncertainty. Although not predictable, if the variations of slopes within spatial scale of soundings are large ($>$±5°), the angular shape of specular region of the backscatter will be distorted (Figure 16). For a given beam, the incidence angle is provided by the best available bathymetry. However, the variations in the local slope can introduce uncertainty in the incidence angle and therefore, while binning backscatter values, the backscatter values from other incidence angles will be averaged together, introducing the uncertainty shown in Figure 16. Even if the highest possible resolution of bathymetry is used to compute the seafloor slope, there will still remain some level of unresolved seafloor slope. However, it is concluded that based on the requirement of taking averages over an angular bin, these unresolved seafloor slopes are problematic only in some cases of the specular region of the backscatter. Binning backscatter data to a resolution similar to the bathymetry grid offers a solution to reduce the impacts of variations in local seafloor slopes. Also, this can provide a practical guideline for applications where selecting a spatial scale of the backscatter strength is required. For example, Buscombe et al. [32] recently suggested spectral filtering to remove high frequency noise from the backscatter data. The low-pass filtered backscatter data was suggested to be representative of the underlying sediment. The selection of spectral filtering parameters can, however, be subjective. Using the bathymetry resolution scale as a guide to the binning size can provide a quantitative and practical means to select a low-pass filter for the backscatter data.

4.1.3. Impact of Bathymetric Uncertainty

Bathymetric uncertainty was shown to cause moderate uncertainty (<3°) in the seafloor slope (Figure 14). This was based on the empirical uncertainty derived for one case-study from CUBE and using a Monte Carlo simulation to run 50 iterations of slope estimation for the bathymetric grid at a 1 m grid cell size. The bathymetric uncertainty, however, will vary for different depths as well as within the swath extent. Multibeam surveys are usually run to comply with bathymetric uncertainty guidelines from the International Hydrographic Office (IHO) [33]; however, modern multibeam systems often surpass these guidelines by far. The total vertical uncertainty (TVU) of soundings is required to be better than the IHO Special Order for shallow water surveys (<40 m) and better than the IHO Order 1 (<100 m) or IHO Order 2 (>100 m) for deep water surveys. The vertical uncertainty, ${\sigma }_z$ (at 95% confidence), is depth-dependent and defined by IHO as [33]:

$${\sigma }_z=\pm \sqrt{a^2+{(b\times d)}^2}$$

(17)

where a is the portion of uncertainty that does not vary with depth, and b is a coefficient that represents the portion of uncertainty that varies with depth, d. For local slope calculations, the soundings are not expected to be affected by non-random (stable) errors (such as biases caused by the system installation or by the sound velocity profile), therefore, only random relative vertical uncertainties should be considered. Values of a and b for the IHO Special Order, Order 2, and Order 1 are provided in

Table 5. Parameters used to estimate vertical uncertainty using IHO [33] uncertainty guidelines for various depths.

	Water Depth (d)	a	b (% d)	${\mathit{\sigma }}_{\mathit{z}}$
Special Order	20 m	0.25	0.0075 (0.75)	0.29
Order 1	50 m	0.5	0.013 (1.3)	0.82
Order 2	100 m	1	0.023 (2.3)	2.50

and can be used to estimate the worst case of vertical uncertainty.

Using the vertical uncertainties provided in Table 5 to perturb the soundings, Monte Carlo iterations applied to Equation (10) can provide slope perturbation estimates as a function of the horizontal spacing of the soundings. The resulting standard deviation of the slope estimates shows that bathymetric data (even complying with IHO Special Order standards) provide very high standard deviations in the estimated slopes (Figure 17). For example, for 0.5 m sounding spacing, the Special Order survey shows a local slope uncertainty exceeding 35°. The dependence of the slope uncertainty on the sounding spacing, as observed in Figure 17, is due to the fact that the same vertical uncertainty will result in higher errors in the computed slope from smaller sounding spacings due to the reduction in the denominator in the slope estimation, from Equation (10). However, for modern MBES, the depth uncertainty is well controlled, and several studies have suggested that the vertical random uncertainty will surpass the Special Order requirements. For example, Marks and Smith [34] found that 95% of multibeam sonar surveys archived at the US National Oceanic and Atmospheric Administration (NOAA) national archives show a repeatable depth within 0.47% of the water depth (to be compared to the Special Order parameter, b = 0.75%). Similarly, performance testing of several new systems (e.g., [35,36]) have shown that depth uncertainty is constrained well below the IHO standards. Estimates of slope uncertainties with Total Vertical Uncertainty (TVUs) of 0.1% and 0.47% of the water depth are plotted in Figure 17 for comparison and show that slope uncertainty can be reduced to manageable levels if vertical uncertainty is constrained well.

The TVU for multibeam soundings strongly depends on the angle with respect to the vertical (tilt angle). For the demonstration, one may adopt a model to approximate the depth dependent TVU of the form:

$$\frac{\delta z}{d}=\sqrt{{\alpha }^2+{\beta }^2{\tan }^2\theta }.$$

(18)

Using typical values of 0.001 for $\alpha$ and 0.003 for $\beta$ provide TVU estimates that replicate approximately the uncertainty observed with respect to the tilt angle [35] (0.1% d near nadir to ~1% d at outer beams). To estimate slope uncertainty using the above TVU estimate, the across-track sounding spacing is assumed based on equidistant beams spread over the angular swath. The across-track distance between the soundings as a function of depth (d) can then be approximated as:

$$\Delta x/d=2\tan {\theta }_M/N_b$$

(19)

where ${\theta }_M$ is the maximum tilt angle and $N_b$ is the total number of beams. Using ${\theta }_M$ as 75° and $N_b$ as 160, the depth-dependent sounding spacing can be computed (19). Considering the vertical uncertainty (18) and sounding spacing (19) are both linearly dependent on the depth, d, the slope uncertainty is finally depth-independent and can be estimated by computing the standard deviation of Monte-Carlo iterations of the slope computation (10). The standard deviation of the slope using the likely TVU and beam spacing (with a 0.1 m position uncertainty) shows the resulting dependence on the tilt angle (Figure 18). For this configuration, the slope uncertainty varies from 3° at nadir to 8° at the swath end. Note that the slope uncertainty can be further reduced by computing slope over larger across-track intervals, at the expense of spatial resolution (Section 3.3); Figure 18 illustrates the reduction in slope uncertainty when increasing Δx by a factor 2 and 4. In conclusion, bathymetric surveys strictly complying with IHO standards can still result in prohibitive uncertainty in local seafloor slopes usable for backscatter computation; however, modern MBES can outperform IHO bathymetric standards by many orders and having an optimal TVU will enable slope to be computed with high (<6°) uncertainty that can be further reduced to moderate levels (<3° as defined in Table 3) if the slope is computed over multiple soundings.

4.1.4. Slope Uncertainty vs. Grid Resolution and Computation Approach

The choice of the seafloor slope computation algorithm does not contribute significantly to the seafloor slope uncertainty. Two computation algorithms were tested during this study. The differences between the results of these two algorithms were observed to be small (<1°). Other studies have assessed the impacts of the choice of slope algorithms and have determined that this effect is not as pronounced as the choice of the spatial scale [37,38]. In conclusion, considering that robust computational algorithms are available, the choice of computation algorithm has limited impact on the seafloor slope uncertainty.

The practical question of assessing uncertainty of the seafloor slope by end-users for seafloor backscatter correction still requires the development of tools accounting for the computation algorithm, scale, and bathymetric uncertainty. An important realization is that the standard deviation of the slope (computed at a single scale) does not represent the actual practical uncertainty in the slope estimation. For example, while considering smaller grid sizes (or small sounding spacing as in the case of beam bathymetry), large values of slope uncertainty are realized. For using coarser grid cell sizes (or large distances between beam soundings), the computed standard deviation of the slope is reduced significantly. This falsely implies that coarser grid resolution results in lower slope uncertainty. An alternate approach to quantify uncertainty is to estimate the seafloor slope at various (possible) resolutions then compare these to estimate the uncertainty of the insonified area and incidence angle. The majority of end-users rely on commercial software packages to process MBES backscatter data [4]. Therefore, they have limited choices of which method to choose and too often these choices are not transparent. In commercially-available backscatter processing tools, no uncertainty in the seafloor slope (either from sounding uncertainty or from the choice of the grid-cell size and methodology) is explicitly defined in the final results. Therefore, validation of backscatter results without the availability of their uncertainties becomes a challenge for end-users. Software developers are hence urged to provide more details about their data processing algorithms and ideally to incorporate uncertainty estimators. These details can be captured effectively in the form of meta data, which should include among the basic information about the MBES (e.g., frequency), the processing parameters and methods used for corrections. Comparing and validating processed backscatter results and providing detailed data processing steps (including insonified area and incidence angles) have been recommended by the ad hoc Backscatter Working Group [18]; this approach will help end-users in comparing and contrasting the effects of the processing options, and hence in optimizing the methodology and the spatial scale of slope correction.

5. Conclusions

Uncertainty quantification is a complex and important part of seafloor acoustic remote sensing that is integral to the repeatability of the processing and interpretation of backscatter data. The issues highlighted in this study relate to the computation of seafloor slopes and their impact on the incident angle and insonified area calculations. The magnitude of seafloor slope uncertainty is impacted by the uncertainty in the measurement of seafloor topography, the methods used to model the topography, and the spatial scale of the bathymetry used to compute seafloor slopes. The order of the magnitude of uncertainty expected from each source was identified (Table 4), showing that the flat-seafloor assumption often used during data acquisition is justified only by the real-time computation constraints and should be avoided in post-processing. Fortunately, software processing tools now enable end-users to correct for the seafloor slope. However, the corrections still can suffer from large uncertainties if the highest possible resolution of the seafloor slope is not used. The spatial scale of the bathymetric data dictates the scale at which backscatter data can be accurately estimated; consequently, the backscatter values should be averaged at the scale of the bathymetry used for slope estimation. It is hoped that by making the issues of the slope, incidence angle, and insonified area more explicit in post-processing operations, future studies of seafloor backscatter will incorporate uncertainty in their analyses of seafloor backscatter strength, end-users will pay more attention to these issues during data processing and interpretation, and software developers will provide processing tools with more accurate compensations of the slope effects.