Modeling the effect of random roughness on synthetic aperture sonar image statistics

Synthetic aperture sonar (SAS) is rapidly becoming a standard tool for seafloor imaging (Hagemann, 1980; Hansen, 2011), seafloor characterization (Williams, 2015; Zare et al., 2017), and target detection (Sternlicht et al., 2016; Williams, 2018). An understanding of the physical processes affecting the acoustic scattering statistics of SAS images is a vital step toward fully utilizing the data produced by these systems for remote sensing, seafloor segmentation, automated target detection, or sonar performance prediction applications. In regions with homogeneous geo-acoustic properties, the local texture seen in images from side-looking systems, such as sidescan sonar or SAS, are indicative of local, small-scale, seafloor slope variations. An example SAS image of a randomly rough sandy seafloor taken in 2015 off the coast of northwest, Florida, using the Naval Surface Warfare Centre's (NSWC) experimental SAS system (Sternlicht et al., 2016) is shown in Fig. 1 and shows increasing variation in intensity as the range from the sonar increases. These increasingly strong intensity fluctuations as a function of range are due to the larger variation in scattering strength versus angle that exists at lower mean grazing angles. These fluctuations in intensity will modulate the imaging speckle, strongly influencing the overall statistical characteristics of SAS images (Lyons et al., 2010).

The continuous variation in scattering strength produced by a random slope field can be treated as an intensity scaling on the image speckle that is produced by the coherent SAS imaging process. Speckle is produced by the coherent summation of scatterers within the system's resolution cell, defined by the ambiguity function of the system (consisting of the pulse, matched filter, and array processing) (Abraham, 2019). If the scatterers can be assumed to be random, and no coherent targets or target-like bottom features (rather than diffuse reverberation) are present, then the process can be assumed to be a multiplicative random process, with the power of the scattered field scaling the speckle component (Oliver and Quegan, 2004). The overall statistics at any position in a SAS image can be expressed as the product of speckle noise and the scaling due to the underlying seafloor scattering cross section:

In Eq. (1), Z(r, x) represents the speckle intensity field for each pixel location in an image and a(r, x) represents a modulating process that captures the effect of the random roughness-induced intensity variation; r and x, respectively, represent the down-range (or simply range) and cross-range (or along-track) image dimensions. In Eq. (1), Y(r, x) is the matched-filtered and beam-formed intensity, so that speckle following a Rayleigh-distributed envelope would produce an exponentially distributed intensity (Lyons and Abraham, 1999). As previously noted, the scaling function a(r, x) is the acoustic expression of scattering from areas larger than the system's spatial resolution and is a function of the seafloor slope field. This scaling function is, effectively, a modulation of the scattering cross section, $σs(θ)$⁠, caused by local grazing angle changes due to variations of seafloor slope, which is treated as a random process. In previous work, a(r, x) was treated as a deterministic process caused by changes in slope due to ripples (Lyons et al., 2010). The cross section term can be calculated via empirical or approximate models of seafloor interface scattering, such as Lambert's law (Lambert, 1760), perturbation theory (Kuo, 1964), or the small-slope approximation (Voronovich, 1985).

In Sec. II, we present a model based on the composite roughness approximation [a combination of the Kirchhoff approximation and perturbation theory (McDaniel and Gorman, 1983)] which has been developed to predict the effect of random, power-law, roughness on the overall SAS image statistics. Changes in image statistics caused by roughness are quantified in terms of the relative intensity variance, or scintillation index (SI). Factors influencing the SI include the slope variance, geo-acoustic properties of the seafloor, such as the ratio of sediment sound speed to bottom water sound speed, the probability density function describing the speckle, and the signal-to-noise ratio. Note that we are treating the scattered field from the seafloor as signal, and any additive components that are incoherent with the signal are treated as noise, which includes ambient noise and multipath contributions (Cook and Brown, 2018; Johnson et al., 2009). Example model-data comparisons will be shown for three different SAS systems from three different sites: (1) SAS images taken in 2010 off the coast of Tellaro, Italy, by the NATO Undersea Research Centre, La Spezia, Italy (now the NATO Centre for Maritime Research and Experimentation) using the 300 kHz Minehunting Unmanned underwater vehicle for Shallow water Covert Littoral Expeditions (MUSCLE). SAS system (Bellettini and Pinto, 2009); (2) data collected in 2013 off the northwest coast of Elba Island, Italy, with the Norwegian Defence Research Establishment's (FFI) 100 kHz Hi resolution Interferometric Synthetic Aperture Sonar (HISAS) system (Fossum et al., 2008) mounted on a High-Precision Untethered Geosurvey and Inspection system (HUGIN) autonomous underwater vehicle (AUV) (Hagen et al., 1999); and (3) data collected in 2015 off the coast of northwest Florida, using the NSWC's experimental SAS system (Sternlicht et al., 2016). Comparisons between parameter estimates obtained from high-resolution SAS data collected in these experiments and historical ground truth will be also used to illustrate the potential use of the model for estimating roughness parameters, such as root mean square slope or height.

For simplicity, we consider only the range dimension when relating the seafloor slope field to image statistics. This approximation is accurate to first-order in root mean square (rms) slope, and this simplification was also performed in previous literature on the composite roughness approximation (Jackson et al., 1986). The geometry for our problem is shown in Fig. 2, defining the true incident grazing angle interrogated by the sonar system, θ, the mean grazing angle, θ₀, and the local slope angle, $ϕ$⁠. The scattering strength will vary about its mean value as a function of range based on the relative (or local) grazing angle at a given range. The concept of larger-scale slopes modulating the scattering from smaller-scale roughness is the foundation of the composite roughness theory for seafloor backscatter (Jackson et al., 1986; McDaniel and Gorman, 1983), that has been used to predict the scattering cross section. Here, we use the composite roughness model to compute higher-order moments of the intensity. Based on Eq. (1), the intensity function at any position in the image, Y(r, x), is the product of speckle intensity and the scattering cross section evaluated at the true grazing angle with respect to the seafloor, $θ=θ0+ϕ$⁠, at range, r:

where $σs(θ)$ is the scattering cross section of the small-scale roughness (within a resolution cell), $ϕ$ is the seafloor slope angle with zero mean, and θ₀ is the nominal (flat-interface) grazing angle at the range of interest. The small roughness perturbation method is used for σ_s, and formulas, approximations, ranges of validity for this model can be found in Jackson and Richardson (2007). In general, this model should be valid for grazing angles and frequencies appropriate for synthetic aperture sonar systems.

We quantify the changes in image statistics versus range caused by random roughness in terms of the normalized intensity variance, or scintillation index (SI) (Jackson and Richardson, 2007). Higher values of SI are indicative of heavier-tailed scattered amplitude distributions and a value of 1 signifies a Rayleigh distribution. The scintillation index,

representative of the image statistics at a given range, can be calculated in our case via knowledge of the intensity moments, m₁ and m₂, of Y(r). Additive noise will also exist for any realistic sonar system. As derived in  Appendix A, the SI of a measurement containing signal in the presence of noise can be described as a modified form of Eq. (3):

This formulation depends on the range-dependent noise-to-signal power ratio, $δ(r)$⁠, which is equal to

The noise-to-signal power ratio is parameterized by the equivalent noise angle, θ_n, the angle at which the noise term and the seafloor scattering cross section are equal. Other parameters in Eq. (5) are defined in  Appendix A. This parameterization allows us to cancel constant terms in the sonar equation for signal and for noise. Local slope is treated as a continuous variable to yield the intensity moments in terms of formal expectations,

where $fs(ϕ)$ is the slope distribution, assumed to be Gaussian with variance $σg2$⁠:

Equation (6) results from using the high-frequency Kirchhoff approximation on the large-scale roughness in a similar fashion as in the composite roughness approximation. Equation (6) is integrated over all slope (positive and negative). The Gaussian assumption for slope follows from the assumption that the seafloor relief is a two-dimensional Gaussian random process as has been assumed in previous texts (e.g., Jackson and Richardson, 2007). Experimental evidence for the assumption of a Gaussian distributed height field can be found in Berkson and Matthews (2009) and Stanic et al. (1989). The underlying speckle statistics are assumed to follow a K-distribution (Abraham and Lyons, 2002; Jakeman, 1980),

with shape parameter α_K, scale parameter λ, and mean power $αKλ$ (which we normalize to 1 in our analysis). In Eq. (8), K is the Basset or MacDonald function (i.e., a modified Bessel function of the second kind) (Oldham et al., 2009) and Γ is the gamma function. Note that Eq. (8) describes K-distributed intensity values which differs from the form found in other articles referenced in this paper (Abraham and Lyons, 2002; Lyons and Abraham, 1999; Lyons et al., 2009; Lyons et al., 2010), which describe K-distributed envelope values. Converting from intensity to envelope would result in the same functional form as given in other references. The shape parameter allows the distribution to provide good fit to a wide range of data. As α_K tends to infinity, the K-distribution tends to a Rayleigh-distributed envelope.

This framework for modeling intensity statistics is similar in spirit to the procedure outlined in Hellequin et al. (2003) who looked at the effects of random seafloor slope on the angular response of multibeam sonar backscatter statistics. We differ from that study, however, in that we perform a numerical integration of the moments in Eq. (3) via Eq. (6) (which only requires a numerical evaluation of the integral over $ϕ$ as for the similar derivation given in Lyons et al. (2010). We also use perturbation theory (Jackson and Richardson, 2007) instead of making an assumption that scattering strength versus grazing angle follows the empirical Lambert's law (Lambert, 1760). The effects of using the more realistic scattering theory will be examined in Sec. III. Also, of note are previous efforts (Lyons et al., 2016; Lyons et al., 2013) that use this basic framework, but do not include the effect of noise, nor test for suitability of the K distribution for the unmodulated envelope probability density function (pdf). We also note the recent work of Olson and Lyons (2021) which used numerical simulations to study the effect of slope modulation on scattering from one-dimensional surfaces. Using the integral equation technique, which is essentially an exact result, Olson and Lyons (2021) demonstrated that intensity fluctuations from a stationary rough surface with Gaussian statistics were due in large part to slope modulations.

We do not include in our modeling, but acknowledge that for very small grazing angles, multiple scattering and shadowing can become important (Liszka and McCoy, 1982). Multiple scattering may become important when the grazing angles are comparable to twice the rms slope and for shadowing when the grazing angle is comparable to the rms slope (Thorsos, 1988). Additionally, volume scattering is often an important contributor to seafloor scattering (Jackson and Richardson, 2007), but is ignored here since our focus in this work is for frequencies on the order of 100 kHz and low grazing angles. For acoustically softer sediments, such as silts and clays (i.e., those without a critical angle), volume scattering may be an important factor to include in any modeling as it may diminish the effects described in this paper depending on the level of volume scatter.

From the previous equations in this section, if a(r) is constant, then the measured SI is that of the underlying distribution. The underlying speckle may itself be non-Rayleigh as noted in Lyons et al. (2010), which is observed in the images that appear in this article in regions where there is no expected modulation due to the random slope fields (moderate grazing angles where scattering strength is fairly constant as a function of angle). If a(r) is not constant, then the image-level statistics will have a scintillation index greater than that of the underlying distribution. At long ranges (corresponding to small grazing angles), the signal produced by scattering from the seafloor may be small compared with additive noise, and SI will then be dominated by the $δ(r)$ terms. For this case of Gaussian additive noise, the SI asymptotically approaches unity. It should also be noted that when comparing these results to real data, system calibration is not necessary because any system dependent parameters, such as vertical beam pattern, would appear in both the numerator and denominator of Eq. (3) for a given range.

Having established a theoretical framework within which to interpret the effects of roughness-induced intensity scaling of the underlying speckle on the statistics of SAS images, we next investigate the applicability of these results by comparison to three datasets. Our first comparison is to an image collected with the 100 kHz HISAS system (Saebo et al., 2007), taken off of Elba Island, Italy, in 2013. The second comparison is to an image collected with the 300 kHz MUSCLE SAS (Bellettini and Pinto, 2009), taken off the coast of Tellaro, Italy, in 2010. The third comparison is to an image collected using experimental SAS system (Sternlicht et al., 2016) obtained off the northwest coast of Florida. The SI estimated from each image is compared to predictions based on the scaled-speckle model developed in Sec. II, namely Eq. (4). The MUSCLE SAS transmitted 60 kHz bandwidth signals at a center frequency of 300 kHz and the HISAS system transmitted 30 kHz bandwidth signals at a center frequency of 100 kHz. Images formed from the collected MUSCLE data had a resolution of approximately 1.5 cm in range × 2.5 cm along-track and images formed with HISAS data had a resolution of approximately 3.5 cm × 3.5 cm. The NSWC system produced data of similar resolution (i.e., on the order of cm) to both the HISAS and MUSCLE systems.

The data from all three of the sites used for model-data comparisons were obtained on uniform and homogeneous seafloor areas of fine sands with varying degrees of topographic roughness, with the Elba experiment conducted in approximately 42 m water depth, the Tellaro experiment in 17 m water depth, and the Florida experiment in approximately 16 m water depth. Although multipath, if present, would impact estimates of SI, multipath was not an issue for the datasets examined as part of this work. Additionally, multipath would mimic additive noise causing estimates if SI to tend back toward smaller values, an effect that we capture via Eq. (4) and Eq. (5). To isolate the effect of roughness changes in our model-data comparisons, SAS images were chosen with little or no variability within an image due to changes in sediment geo-acoustic properties. Image selection was accomplished by visual inspection of the image at moderate grazing angles where the scattering cross section is insensitive to small changes in grazing angle (as discussed in  Appendix B). Areas very close to the specular direction were excluded, due to the rapid variation in scattering strength with angle, as well as beam pattern sidelobe artifacts in that region.

Geo-acoustic input parameters for the perturbation theory-based model for the same general vicinity as the Elba Island HISAS data collection site were published in Maguer et al. (2000) and Lyons et al. (2009). Inputs to the scattering model used for comparisons with the MUSCLE SAS data are from Pouliquen and Lyons (2002) and were obtained close to the location off of Tellaro, Italy, where the SAS image data were obtained, and at the same water depth. Inputs to the model used for the northwest Florida dataset were from Williams et al. (2009). The perturbation-theory scattering model used in the SI model requires a roughness power spectrum as input, and we use the power-law form,

which has been found to provide a good fit to a wide range of seafloor roughness measurements (Jackson and Richardson, 2007). In this equation, $k$ is a two-dimensional wave vector with magnitude $k=|k|$⁠, β is the spectral strength parameter, and γ is the spectral slope parameter. The seafloor slope variance over the scale of a SAS image used in Eq. (7) can be obtained by integrating the power-law roughness slope spectrum between an outer scale k_L, and an inner scale k_c, as outlined in Jackson et al. (1986) 

where $kc=2π/dx$ is inversely related to the smallest scale in the image, dx (i.e., dx is the sonar resolution or pixel size if the image is not oversampled), and k_L is inversely related to the largest scale in the image (i.e., the total along-track distance of an image as we will calculate SI on individual range lines). Note that Jackson et al. (1986) uses a different form of the slope pdf than our Eq. (7). Assuming isotropy yields

Note that the rms slope depends both on the sonar resolution and the image size, although it is less sensitive to image size, since k_c is typically much larger than k_L. Shape parameters are estimated from the data as outlined in  Appendix B.

For the comparisons that follow, we hold the spectral slope, γ, constant and vary the spectral strength, β, as a free parameter to fit SI versus range. The noise parameter, θ_n is also fit to the data. The dependence of SI on range allows us to fit the rms slope (via) in an unambiguous fashion via the following steps:

• Estimate the free parameters in the K-distribution from the near range image part that has relatively larger grazing angles and is the most unaffected by slope modulation.

• Estimate β (and therefore σ_g) by matching the theoretical model to the data, Eq. (3) ignoring noise, using moderate ranges where $SI$ increases. Eq. (6) is numerically solved during this fit.

• The additive noise parameter, θ_n, is estimated by matching the long-range data, where SI tends to flatten with range.

Table I gives parameters used in estimating SI for each of the study sites. In this table, ρ is the sediment-water density ratio, ν is the sediment-water sound speed ratio, and η is the loss parameter of the sediment, defined as the ratio of the imaginary to real components of the wavenumber in the sediment. The complex sound speed ratio is $ν/(1+iη)$⁠. More detailed descriptions of these geo-acoustic parameters can be found in Jackson and Richardson (2007).

Figure 3 shows a SAS image and estimates of SI versus range for the HISAS Elba Island dataset. Each SI estimate in range was formed for a single range line only, i.e., there is no sharing of data across range cells. The experimental SI versus range estimates have had outliers removed and have been smoothed with a 21-point averaging window. The outlier filter simply compares each data point to a local average and removes the data point (without replacement) if the difference between the two is greater than a given threshold (a threshold of 1 was used in our analysis). The SI versus range estimates with outliers removed and then smoothed were then used to fit the model presented in Sec. II for model-data comparisons. Given that 21 points in range represent only a small change in angle, smoothing is expected to have little effect on model fits. The SI is seen to increase dramatically as a function of range (i.e., increase as the mean grazing angle decreases from approximately 22° at 40 m range to approximately 6° at 150 m range). Although it is the non-Rayleigh speckle statistics that cause the value nearer nadir to be larger than 1, it should be noted that it is not the underlying speckle statistics that are the dominant cause in the overall increase in the SI versus range, but the increasing slope of the scattering cross section versus grazing angle at further ranges (or smaller grazing angles). We also note here that with the implicit assumption that the seafloor height field is isotropic and statistically stationary, there will be no dependence of SI on the azimuthal look direction of the sonar system.

An illustration of this effect for various seafloor classes is shown in Fig. 4. The top figure shows scattering strength as a function of grazing angle for various sediment types, whereas the bottom plot shows the absolute value of the gradient of the backscattering strength. The input parameters to this model were based on recommendations from Jackson (2000) and Richardson and Jackson (2017). Volume scattering is included in each of the curves but contributes insignificantly to the scattering strength curve, except for silt. To compare with our work, the medium sand parameters from this figure are closest to the geoacoustic parameters used for the model-data comparisons presented in Figs. 3 and 5–8. Taking the gradient of the logarithm of the scattering cross section results in an estimate of the gradient of the cross section results divided by the cross section, which is the important parameter for this work. Different seafloor types have relatively consistent, small gradients at moderate grazing angles, and all increase dramatically at small angles.

The scaled-speckle model developed in the Sec. II was used to predict the SI using the rms seafloor slope as a free parameter and was fit to the experimental estimates of SI using known system geometries. A K-distribution shape parameter value of 4.0 and equivalent noise angle of 0 was used in the model calculation for this comparison. The values of SI at the smallest ranges in Fig. 3 (i.e., with the largest mean grazing angles), where scattering strength is close to constant as a function of angle, are almost solely produced by the speckle statistics. An rms slope value of 4.9° was found to provide the best fit to the experimental data displayed in Fig. 3.

To illustrate the effect of the choice of scattering model, SI was also calculated with the empirical Lambertian scattering model. The model results displayed in Fig. 3 show that SI increases much more quickly as range increases (as grazing angle decreases) for predictions based on perturbation theory than for the Lambertian model. The intensity variability seen in an image at a specific range is directly related to the shape of the scattering strength curve, i.e., the rate of change (and higher-order derivatives) of scattering strength as a function of angle. The approximately $sin4$ dependence of scattering strength versus grazing angle at low angles for perturbation theory yields a larger SI for smaller grazing angles than the Lambertian case (with the same rms slope) while the relatively constant scattering strength versus angle near the critical angle of 27° yields a flatter curve of lower SI equal to that of the underlying speckle.

The dominant factors affecting the SI versus range are the K-distribution shape parameter of the speckle and the rms seafloor slope (which we parameterize using the spectral strength). Figures 5 and 6 display the impact of these two parameters on SI. In these figures, comparisons are made with the same HISAS-derived SI estimate versus range as displayed in Fig. 3. In Fig. 5, the shape parameter is varied and rms slope held constant. It can be seen in Fig. 5 that the influence of the shape parameter on SI asymptotes relatively quickly to having little impact as shape parameter increases. In Fig. 6, the shape parameter is held constant, while rms slope is varied. Changes in shape parameter for the most part only shift the SI versus range curve in its entirety to higher or lower values, while the spectral strength (or, equivalently rms seafloor slope) changes the slope of the SI versus range curve, which appears to reach an asymptote. The dependence of the SI curve at close ranges on the shape parameter (as discussed in  Appendix B) is evident when considering Figs. 5 and 6.

The second dataset used to explore the effects of random roughness on SAS image statistics consisted of data taken in 2010 off the coast of Tellaro, Italy. An example image is shown in Fig. 7. The SAS image in Fig. 7 clearly displays more large-scale variation in intensity from about 90 m in range onwards than the images in Fig. 1 or Fig. 3, which is caused in this case by larger undulations in seafloor slope. These very-large-scale topographic undulations (and associated intensity variation) are not accounted for in our model. The bottom plot in Fig. 7 shows estimates of SI versus range estimated from the MUSCLE dataset that were used to form the image shown in Fig. 7. As for the Elba Island HISAS data, the SI is seen to increase as a function of range away from the sonar (i.e., SI increases as the mean grazing angle decreases from approximately 14° at 40 m range to approximately 4° at 145 m range). The large fluctuations in the SI seen in this figure are due to tilting of the large-scale slope toward and away (along with some subsequent shadowing) from the sonar. It should be noted here that large-scale intensity fluctuations in SAS imagery may also be caused by refractive effects caused by internal wave related features as recently noted by (Hansen et al., 2015). The effects on the SAS image statistics caused by these types of intensity fluctuations would be indistinguishable from those caused by large-scale topographic slope variation.

One aspect that appears in the measured SI for the Tellaro data that was not seen in the HISAS data shown in Fig. 3 is a flattening of the curve at larger ranges. This effect is caused by a loss of signal-to-noise ratio at far ranges (or in areas dominated by shadows) which drives the SI toward lower values (i.e., in the direction of the SI for a Rayleigh distribution) as Gaussian additive noise begins to have an impact on the sample distribution. The noise angle, defined as the angle at which the noise term and the seafloor scattering cross section are equal, was set to 2.8° for the comparison shown in Fig. 7 (an equivalent range of approximately 200 m).

For our third example dataset, we show results for data taken with the NSWC system off the northwest coast of Florida. Measured and modeled SI for this site (and system) are shown in Fig. 8. In this example, both port and starboard data and model results are given and agree with the previous two datasets in showing a rapid increase in SI versus range. The two-sided results shown in this figure highlight another interesting aspect of roughness-induced SI versus range; the asymmetry caused by a tilt angle of the sonar collecting the SAS data with respect to the seafloor. This tilt can be caused by either a tilt in the vehicle housing the port and starboard sonar systems or a slight image-wide slope in the seafloor. For the comparison shown in Fig. 8, we have included a 0.5° relative tilt between the sonars and the seafloor in the model.

Roughness values are often reported in terms of rms roughness (over a specified length or area). These values can be compared to the values resulting from our model-data fitting to lend confidence in our methods. Statistical properties of the seafloor roughness, such as rms height and slope, can be related if a model of seafloor roughness is assumed. If a seafloor exhibits power-law roughness spectra (as the model developed here does), rms height, h, can be obtained from rms slope as outlined in Jackson et al. (1986),

Under the assumption of power-law roughness, the rms seafloor height was found to be to 7 cm for the Elba Island, Italy, dataset, 5 cm for the Tellaro, Italy, dataset, and 2.4 cm for the northwest Florida dataset. For self-affine surfaces, such as those described by a power law, surface measures, such as rms height or roughness, h, will depend on largest observation length, l, and obey a scaling relationship of the form (Summers et al., 2007)

Jackson and Richardson (2007) presented stereo-photogrammetry–based measurements of rms roughness of 0.2–1.0 cm for fine and medium sands measured over scales on the order of 1 m. Using Eq. (7) with a spectral exponent of 3 (approximate value for fine and medium sands is based on Table VI.1 in Jackson and Richardson (2007)], these historical values of h would yield rms roughness values measurements of approximately 2.5–12 cm over the much larger 50–150 m scale of the SAS images, values which bracket the 2.5–7 cm rms roughness predicted by the scaled-speckle model for the three sites presented in this study. Seafloor roughness data measured by divers or estimated via stereo photogrammetry methods at Elba Island, Italy (Maguer et al., 2000), Tellaro, Italy (Pouliquen and Lyons, 2002), and northwest Florida (Williams et al., 2009), all close to the sites of the data presented in this paper, reported rms roughness on the order of what our results suggest via fitting of the SI model (i.e., on the order of 3–10 cm). Although roughness will not be temporally constant, the mechanisms generating roughness at a given location (and depth) with given sediment properties (e.g., grain size) may be relatively stable. The relative agreement between reported roughness values and our own estimates suggests that the parameters that best fit the data presented here are reasonable. This model could potentially be used over larger areas to invert for seafloor roughness properties, which are useful for hydrodynamic modeling, sediment transport, and acoustic reverberation–an important component of sonar performance models.

In this work, we have presented a model to predict the impact of random seafloor roughness on SAS image statistics. In the model that was developed, the variations in scattering strength produced by changes in seafloor slope were treated as a scaling of the coherent imaging speckle produced by a synthetic aperture sonar. The changes in image statistics were quantified in terms of a common complexity metric, the SI. For the three experimental datasets examined, roughness caused a dramatic effect on the statistics of the images, with increasingly larger SI (i.e., heavier-tailed distributions) as range away from the sonar increased. SI estimates from SAS data showed good agreement with model predictions, showing sensitivity to the roughness spectral strength (i.e., the mean square slope). This sensitivity could allow seafloor roughness parameters, such as rms roughness, to be inverted from SAS image data as in Chen et al. (2015). The range dependence of the effects of roughness on image statistics should also be considered when setting thresholds in automatic target recognition (ATR) detectors and when using scintillation-related metrics for seafloor characterization or image segmentation (e.g., those which use lacunarity (Williams, 2015), which is equivalent to scintillation index when using intensity values). Last, as the effect of the seafloor slope distribution on image statistics should increase with range, a trend back toward a value of 1, or even a flattening of the SI curve (as seen in Fig. 7) will be an indication of a lessening of image quality as noise effects begin to affect the image.

The authors gratefully acknowledge Warren Fox, Hans Groen, and others at the NATO Undersea Research Centre [now the Centre for Maritime Research and Experimentation (CMRE)], La Spezia, Italy, for allowing data from the AMiCa'10 sea trial to be used for this study. The authors also thank the CMRE for organizing the MANEX'13 trials onboard NRV Alliance and the Norwegian Defence Research Establishment for gathering the HISAS data used here with a HUGIN AUV during the trials. The authors thank Darrell Jackson from Applied Physics Laboratory at University of Washington, Seattle, for providing code to produce the sediment backscatter model used in Fig. 4. Thanks also to the Naval Surface Warfare Center Panama City Division for providing the experimental SAS data that was presented in this paper. A. L. was supported under ONR Grant Nos. N00014-10-1–0051, N00014-10-1–0047, N00014-10-1–0151, and N00014-16-1–2268 and D. O. was supported under ONR Grant No. N00014-16-1–2335.

To deal with the effect of noise, we consider the sum of two independent complex random variables U and V. We treat the seafloor return at any range, r, as V, and the additive noise component as U, and assume its real and imaginary parts are independent Gaussian random variables (i.e., are complex-Gaussian distributed). The SI of the sum of these two random variables can be calculated using Eqs. (6) and (8) from Lyons et al. (2009) 

where $m1X=E[|X|2]$ is the first moment of the intensity (second moment of pressure) of random variable X, and $m2X=E[|X|4]$ is the second moment of intensity (fourth moment of pressure). If we assume that U has Gaussian real and imaginary components, then $m2U=2m1U2$ (Lyons et al., 2009). We introduce the noise-to-signal power ratio, $δ=m1U/m1V$⁠, and are therefore able to express all the U moments in terms of moment of V. Performing these substitutions results in the modified SI

which can be rewritten in terms of the SI without noise, i.e., Eq. (3) as

As a check on these results, we examine three limits. First, the limit of $δ→0$ is when noise is completely absent. This results in the original formula for the SI of V, which it should. Second, we examine $δ→∞$⁠, where the signal is only Gaussian noise. This results in a SI of unity, which is expected. Third, we examine the case where V has Gaussian real and imaginary components. In this case, the SI is unity independently of δ, which is what is expected when two complex circular Gaussian random variables are summed.

To model the noise to signal power, δ, we use the sonar equation for synthetic aperture sonar received signal intensity (Cook and Brown, 2018; Massonnet and Souyris, 2008; Olson et al., 2016)

where C_s includes all constant terms in the sonar equation for the seafloor signal (including source transmit power and element gains), r is the slant range to a given pixel, σ_s is the scattering cross section, A is the ensonified area, G is the increase in signal power after synthetic aperture array beamforming, and $kw″$ is the imaginary part of the complex wavenumber in water. For SAS, A is constant as a function of rang, (except for a weak $1/ cos θ$ dependence which we ignore at small angles), and G is proportional to r², since the number of elements for a given pixel is proportional to range. The image noise intensity, I_N, is given by $IN=CnInr$⁠, where I_n is the noise intensity at a hydrophone (assumed constant and isotropic), and the factor of r is due to the increase in noise power after SAS processing. The constant C_n takes into account constant factors in the noise sonar equation, such as the synthetic aperture weighting function [e.g., Hamming or other apodization function (Harris, 1978)]. Combining the constant terms into a single constant, C, we obtain

Although δ is dimensionless, it depends on range. A convenient way to capture this range dependence is to define an angle at which δ is unity, that is when the signal and noise power are equal. From the vehicle altitude H, we may also define a range at which this occurs, $rn=H/ sin θn$⁠. Thus, we have

Solving for I_n and inserting into Eq. (A5), we obtain

Thus, we may calculate SI using Eqs. (A3) and (A8).

The scaled-speckle model developed in this work requires estimates of the parameters of the K-distribution for the speckle component without the influence of a modulating slope field. For normalized data, the K-distribution model requires only that the shape parameter be estimated. In SAS images, a region usually exists close to nadir that is minimally influenced by seafloor slopes, and which therefore appears in the images as a region of low contrast or low texture. This low-contrast, low-texture region exists because the shape of the scattering strength versus grazing angle curve is relatively flat for angles higher than about 15° (i.e., areas closer to nadir, but not near the specular direction). This area of low texture at close ranges (corresponding to moderate grazing angles) can be seen in Fig. 1 at a range of approximately 5–15 m. SAS images at angles very close to nadir have greater contrast due to the large gradient in scattering strength with angle near the specular direction (Jackson and Richardson, 2007), as well as beam-pattern artifacts. With the assumption of no impact of the slope field at near ranges in SAS imagery, it is possible to use data from this region to estimate the speckle shape parameter.

We acknowledge that there will be some small influence of the slope field on the shape parameter at near ranges that cannot be removed. In fact, numerical tests in (Olson and Lyons, 2021) for pressure-release surfaces revealed that for very high bandwidth pulses (comparable to the center frequency and bandwidth of HISAS data examined here), SI is greater than 1 for both moderate and small grazing angles. The impact of this influence will be to slightly reduce the rms slope that is used in fitting the model to experimental data.

Before parameter estimation, the data in each image were normalized in the range direction to a mean power of 1 to remove the large-scale effects of spreading loss, absorption, and transmit and receive beam patterns. As neighboring data values showed slight correlation, that data were also decimated by a factor of two by using every other pixel in a range line to obtain independent samples. The shape parameter was then estimated range line by range line for 150 lines of data (which amounted to between approximately 2 and 5 m of the data closest to nadir, but away from visible beam pattern artifacts, depending on the SAS system). Estimation of the shape parameter was carried out using a common statistical estimation method, the method of moments, which uses the first and second sample moments (Joughin et al., 1993; Abraham and Lyons, 2010).

To evaluate the validity of the fit of the K-distribution to the experimental data, non-parametric Kolmogorov–Smirnoff (KS) test statistic p-values (Papoulis and Pillai, 2002) were used to compare theoretical and observed exceedance distribution functions (EDFs), which are defined as one minus the cumulative distribution function (CDF). Sample EDFs for the Tellaro site are shown in Fig. 9 as a function of normalized envelope amplitude. EDFs from individual range lines are shown on the figure (thin gray lines) along with their median (thick black line) and compared with both the Rayleigh and K-distribution model EDFs (thick dashed line and thick solid line, respectively). The K-distribution shape parameter used in the model fit is the median value of shape parameter found for all the EDFs. Qualitatively from Fig. 9, the K-distribution model provides a good fit to the SAS data (and a much better fit than the Rayleigh distribution). The results for the other two sites were similar to those seen in Fig. 9. The top plot in Fig. 10 shows shape parameter estimates for the 150 range lines of data closest to nadir (excluding the specular region) for the three sample sites (and systems), with the median values shown as solid lines.

The KS test statistic, D, is the maximum absolute error between the sample and model EDFs evaluated over the region where the sample EDF is greater than $10/n$ (i.e., where the EDF has greater accuracy (Preston and Abraham, 2015) and is given by

where $FY(Y(i))$ is the model CDF being tested, with $Y(i)$ being the ordered samples of the echo intensity, and the second term within the absolute value signs representing the sample EDF evaluated at the ith ordered sample. The p-value is the probability of observing a value greater than d under the null hypothesis–the data are distributed according to $FY(Y(i))$⁠. The p-value of d is utilized to compare the results of the KS test. The p-value is the probability of observing a value greater than d under the null hypothesis; i.e., if d is the observed KS test statistic and D is the random variable formed by Eq. (B1) using the random variable Y distributed according to $FY(Y(i)$⁠, then the p-value is

When the p-value is small, observing a more extreme value for the KS test statistic is not likely under the null hypothesis, and therefore the data do not follow the hypothesized CDF. Likewise, when the p-value is large, the maximum absolute difference between the model CDF and the sample CDF is small, and it is unlikely that the model distribution does not provide a good fit to the data. The bottom plot in Fig. 10 shows the p-values for the 150 range lines of data closest to nadir for the three sample sites along with the p = 0.05 criteria for acceptance shown as a dashed line. The K-distribution was accepted at each site as the median p-value for each of the sites was approximately 0.75.