Stochastic formalism-based seafloor feature discrimination using multifractality of time-dependent acoustic backscatter

Dual-frequency echo-envelope data acquired using the normal-incidence single-beam echosounder system (SBES) have been examined to study its scale invariant properties. The scaling and multifractality of the SBES echo envelopes (at 33 and 210 kHz) were validated by applying a stochastic-based multifractal analysis technique. The analyses carried out substantiate the hierarchy of multiplicative cascade dynamics in the echo envelopes, demonstrating a first-order multifractal phase transition. The resulting scale invariant parameters ( α, C1, andH ) establish gainful information that can facilitate distinctive delineation of the sediment provinces in the central part of the western continental shelf of India. The universal multifractal parameters among the coarse and fine sediments exhibit subtle difference in α andH , whereas the codimension parameter C1 representing the sparseness of the data varies. The C1 values are well clustered at both the acoustic frequencies, demarcating the coarse and fine sediment provinces. Statistically significant correlations are noticeable between the computed C1 values and the ground truth sediment information. The variations in the multifractal parameters and their behavior with respect to the ground truth sediment information are in good corroboration with the previously estimated sediment geoacoustic inversion results obtained at the same locations.


Introduction
Acoustic remote sensing methods using the normalincidence single-beam (SBES) and multi-beam echosounder system (MBES) are mainly concerned with identifying, classifying and mapping surficial geological features of the seafloor.These methods are well recognized as a useful tool in oceanography to characterize the seafloor over a wide area and facilitate preliminary geological analyses (Anderson et al., 2008).The seafloor characterization and classification methodologies available in the literature using SBES can be traditionally grouped into two categories, namely modelbased techniques and empirical methods.The model-based techniques often utilize physics-based acoustic backscatter models to characterize the seafloor sediments by optimizing the match between the measured and the modeled signals (Sternlicht and de Moustier, 2003a, b;van Walree et al., 2006;De and Chakraborty, 2011;Snellen et al., 2011).The empirical methods however rely on the statistical and artificial neural network-based approaches to correlate the features of echo signals with the sediment type (van Walree et al., 2005;Chakraborty et al., 2004;De and Chakraborty, 2009;Amiri-Simkooei et al., 2011;Madricardo et al., 2012).
The success of the model-based inversion procedure depends on the scattering theory employed in the forward backscatter model and requires detailed understanding of the scattering mechanism.The study of sound interaction with the seafloor and the corresponding inversion modeling pose a challenging task, particularly with the existence of diversity in the benthic habitat of the area (Holliday, 2007).The scattering process of an acoustic wave is influenced by the presence of benthic fauna responsible for modifying the smallscale morphological features and the density fluctuations within the sediment volume (in addition to the hydrodynamic processes).Incorporation of the number density of biological organisms and their collective activities (i.e., burrowing and home building) into the forward backscatter model complicate the inverse modeling even further.The continuous form of seafloor heterogeneity (due to bioturbation, sediment deposition, or hydrodynamic processes) therefore necessitates the development of versatile and robust statistical techniques to determine the seafloor roughness statistics (Jackson and Richardson, 2007).Accordingly, to further improve the seafloor feature discrimination we introduce an empirical method that uses the scaling and multifractality of the dual-frequency SBES echo envelopes at 33 and 210 kHz.
The "stochastic" multifractal formalism followed herein discriminates the SBES echo envelopes with three fundamental parameters, namely, degree of multifractality α, sparseness C 1 , and degree of smoothness H .In the specific framework of stochastic-based "universal multifractals" (Schertzer and Lovejoy, 1987, 1991, 1997;Lovejoy and Schertzer, 1990), the statistics of the underlying cascade process are completely characterized by the aforementioned scale invariant fundamental parameters (Gagnon et al., 2006).The reason for using a multifractal framework is to build on the fact that the layers of the seafloor imprint a fractal signature on the echo signal along with the self-similarity of sediment ripples of various sizes.Moreover, acoustically soft sediments are penetrated more deeply by acoustic signals and produce longer and corrugated echoes than hard sediments, evidencing fractal structures.Therefore, the estimated multifractal parameters of echo envelopes as a measure of complexity and roughness proffer useful information to improve seafloor feature discrimination.
The remaining part of the work is organized as follows.Section 2 describes the study area, data set descriptions and single-beam data processing methodology.A brief theoretical overview of a universal multifractal framework is also recapped in this section.Section 3 describes the interrelationships among the estimated multifractal parameters in terms of ground truth sediment information, followed by concluding remarks described in Sect. 4.

Acoustic and sediment data acquisition
The dual-frequency (33 and 210 kHz) echo data were acquired over substrates ranging from clayey silt to sand in the central part of the western continental shelf of India, using a hull-mounted normal-incidence RESON-NS 420 SBES.The beam width of the echosounder transducer for 33 and 210 kHz is 20 • and 9 • , respectively, with respective pulse lengths of 0.97 and 0.61 ms.The raw analog output on the receiver circuit board was tapped and connected to a PCL 1712L 12-bit A/D converter with a sampling frequency of 1 MHz.The echoes were stored together with the information of the echosounder adjustments and ship position obtained from the GPS system.Figure 1 shows the 20 locations where the single-beam echoes and sediment data were collected.The seafloor depth of the study area varied between 29 and 109 m (Haris et al., 2011).
The sediment data were collected using a Van Veen grab, covering an area of 0.04 m 2 and penetration of 10 cm, following a standard protocol (at all 20 locations where the echo data were acquired).About 20 g of sediment were taken from each grab sample to carry out the textural analyses using a 4.0 cm diameter core tube.The sediment was repeatedly washed in distilled water until all the chloride ions detectable with 4 % silver nitrate were removed.These samples were treated with 10 % sodium hexametaphosphate and kept overnight for dispersion before being subjected to the grain size analyses.The acquired sediment samples were subjected to wet sieving using a 62 µm sieve to separate the sand from the mud fraction.The size distribution of the mud fraction (< 62 µm) was measured with a Malvern laser particle size analyzer (MASTERSIZER, 2000).The size distribution of the sand fraction was determined using a standard dry sieving method, as it was difficult to maintain uniform suspension of sandy material within the laser particle analyzer.The shelf sediments normally contain shelly material, which had to be sieved prior to measurement by laser diffraction.The mean grain size M ϕ = −log 2 U g /U 0 (where U 0 = 1 mm) was then calculated for each of the sediment data locations.

Single-beam data processing
The recorded echo data were converted from binary to ASCII format within a range of −5 to +5 V. Hilbert transformation was employed to obtain the echo envelope from the echo trace at each location.The shape of the echo envelope is generally influenced by various factors including natural variability of the seafloor, transducer heave, and noise due to echosounder instability (Haris et al., 2012).Therefore, several post-processing steps such as visual check, echo alignment and echo averaging were performed to obtain good averaged echo envelopes to succeed the stochastic multifractalbased analyses.
The first step taken in the post-processing was to select good echo envelopes (with minimum distortion) by removing the saturated and clipped echo envelopes through visual inspection.This was achieved by careful selection of the echo envelopes, characterized by a well-defined initial rise and amplitude followed by a slow decay.In addition, the data having voltage response saturates at ±5 V and those with prominent multi-peaks were discarded.Precisely, an echo envelope displaying only one important peak, with an amplitude between 2.5 and 5 V, was qualified for further processing steps.Due to the transducer heave motion and small variations in seafloor depth over consecutive pings (while recording the backscatter data), initial rise times of the echo envelopes are not the same.Hence, it was imperative to align all echo envelopes before carrying out further processing.The alignment was based on identifying and indexing a temporal feature on an echo envelope.The initial rise time and the time of peak amplitude were considered as an important temporal feature of an echo.After identifying the temporal feature, all echoes within the ensemble were shifted in time to align with the selected feature (Fig. 2).The aligned echo envelopes were averaged to obtain stable acoustic signals to compute the values of universal multifractal parameters associated with different sediment provinces.The echo envelopes were first averaged using 20 successive envelopes with 95 % overlap (in a moving average sense with sequences 1-20, 2-21, and so on, utilizing all the consistent echo envelopes available in the data set).The voltage form of the aligned data was converted to a pressure signal using the hydrophone sensitivity values provided by the transducer manufacturer.The resulting aligned pressure curves were finally ensemble averaged to obtain a representative stable acoustic signal at each location (Fig. 1b).

Depth-dependent correction
Apart from the processing steps described in the preceding subsection, the echo-envelope data require an additional correction for the sonar footprint dimension prior to stochastic multifractal analyses.The footprint diameter enlarges in proportion to the water depth, correspondingly the backscatter area and echo duration increases.As a result, an echo recorded at a greater depth is expanded in time and an echo recorded at a lesser depth is compressed along the time axis (when compared with a reference depth).Consequently, the acoustic returns from the same seafloor sediment type ly-ing at different depths do not have the same shape (De and Chakraborty, 2009).Accordingly, a first-order correction was applied to remove the influence of the depth on the time spread.The time spread of the echo envelope was multiplied by a factor h ref /h, where h ref is a reference depth of 50 m (approximate average of all the spot depths) and h is the depth at the position of the individual echo data (De and Chakraborty, 2010).The procedure followed is equivalent to the depthdependent correction described by van Walree et al. (2005).

Moment scaling function and universal multifractals
Self-similar fractals are scale invariant, i.e., possessing a structure with a basic characteristic of nonscaling.They can be divided into two categories.The first one is the monofractal, having strict geometric self-similarity that can be described with a single fractal dimension.The other is the multifractal that requires a series of fractal spectra rather than a unique fractal dimension.Highly intermittent multifractal fields common in nature are the generic outcome of multiplicative cascade processes dominated by scaling nonlinear interactions.Many geophysical fields have been shown to be multifractal over various ranges (Lovejoy and Schertzer, 2007a).However, in the specific case of echo envelopes, the power-law behavior of the data calls for multifractal analyses to verify the existence of multiscaling.When a multifractal cascade has proceeded over a scale ratio λ = L/ l (L and l representing the largest and smallest timescales in the data), the statistical moments of the conserved multifractal flux (the pressure values of the echo envelope) measured at scale λ follow a power law that can be expressed as (Schertzer andLovejoy, 1987, 1991): where ϕ λ is the scale-by-scale conserved multifractal flux, q is the order of the moment, and K(q) is a nonlinear convex function.K(q) characterizes the scaling of the moments of the ϕ λ , hence it is called the "moment scaling function".With reference to the existence of stable attractive multifrac-tal processes called "universal multifractals" (Schertzer and Lovejoy, 1987, 1991, 1997;Lovejoy and Schertzer, 1990), K(q) can be expressed as: where α and C 1 are the basic parameters characterizing the scaling properties of the multifractal flux ϕ λ .The parameter α is the degree of multifractality and varies from 0 to 2, where α = 0 is the monofractal case and α = 2 is the log normal case.This parameter describes how rapidly the fractal dimensions of the sets at different thresholds vary as they leave the mean singularity.C 1 is the codimension parameter of the set.Low value of C 1 (≈ 0) implies that the field values are close to the mean.C 1 (> 0) indicates that the region making the dominant contribution to the mean is a sparse fractal set such that the vast majority of the field does not contribute (Gagnon et al., 2006).

Fractionally Integrated Flux model
The multiplicative process (the cascade) discussed above generates a scale-by-scale conserved multifractal flux ϕ λ characterized by a moment scaling function K(q).The spectrum of such a conserved flux has an exponent β = 1 − K(2) < 1.In order to discriminate the seafloor echo envelopes (having β ≈ 2), the fractionally integrated flux (FIF) model (Gagnon et al., 2006) has been utilized.The FIF model of the multifractal flux provides the following statistics in relation to the intensity field I λ (pressure values) at scale ratio λ as (Schertzer andLovejoy, 1987, 1991): Here the linear scaling λ −H corresponds to a fractional integral of order H .The parameter H can be designated as a degree of smoothness where a higher H signifies smoother fields.Characterization of seafloor backscattering using the FIF model is difficult to distinguish the underlying cascade dynamics as it involves a convolution due to the exponent H . Therefore resorting to the use of "trace moments" (that directly characterize the conserved multifractal flux ϕ λ ) is necessary so that the differentiation is possible (Gagnon et al., 2006;Chakraborty et al., 2013).The first step to obtain ϕ λ from the intensity field involved the removal of λ −H in Eq. ( 3).This is equivalent to a filtering as in Fourier space with "power law", which is a scale invariant smoothing.On elimination of λ −H , only the underlying conserved multifractal flux ϕ λ is retained.The next step was to examine the scaling of the statistical moments of ϕ λ and compare them with Eq. (1).To this end, we normalized ϕ λ so that the ensemble average of all the samples is < ϕ λ > = 1.Thereafter, the qth power of the samples (pressure values) over the sets (or time interval) of duration l = L/λ was determined.It gives the moments of the normalized multifractal flux for a given value of q.This procedure was performed with different values of q and K(q) was evaluated from the logarithmic slopes (Fig. 3).The multifractality of the intensity field has been validated with nonlinear K(q).From the values of K(q) the parameters C 1 and α were estimated as C 1 = K'(1) and α = K "(1)/C 1 (Stolle et al., 2009;Gires et al., 2013).The values of α and C 1 combined with spectral slope β were utilized to estimate values of H , using the relationship β = 1 + 2 H -K(2).The three universal multifractal parameters (α, C 1 , and H ) computed here determine the statistics of the data at all scales and moments.

Results and discussion
The following sections describe the analyses of the estimated universal multifractal parameters (at 33 and 210 kHz) along with ground truth values of the mean grain size of the seabed sediment.The end results are statistically analyzed and compared with the ground truth data as well as with the previously estimated sediment geoacoustic inversion results (De and Chakraborty, 2011) obtained at the same locations.For simplicity, in the following text, silty-sand and sand sediments will be referred to as coarse sediments (with M ϕ < 4); and clayey-silt and silt sediments will be referred to as fine sediments (with M ϕ > 4).

Sediment distribution
The percentage distribution of sediment composition indicates the presence of four sediment types: clayey silt, silt, silty sand and sand with varied levels of mixing of three textural grades, namely sand, silt, and clay.The important substrate characteristics of the study area have been investigated in detail (Haris et al., 2012), and four distinct sediment provinces were identified from the map generated using the Geographic Information System (GIS) (Fig. 1).The sediment texture was relatively coarse in the deeper depths (60-109 m), whereas fine-grained sediment was found in the shallow depth region (29-54 m).The high percentage distribution of fine-grained sediment in the shallow depth region is governed by a set of sedimentological and hydrodynamic conditions.The study area receives relatively high annual rainfall, which can bring the sediment load through the rivers and discharges to the shallow water regions of the study area.Besides, the shallow region is influenced by an environment with a feeble current that allows fine particles to settle down as compared with the regions of higher depth.These processes might have resulted in the accretion of fine sediment in the shallow depth region.

Multifractal phase transition
Self-organized criticality (SOC) was first introduced by Bak et al. (1987) as an explanation for the 1/f noise detected in various dynamical systems.The SOC is generally evident in multifractal processes along with a multifractal phase transition (Schertzer and Lovejoy, 1992;Schertzer et al., 1993;Hooge et al., 1994;Schmitt et al., 1994;Garrido et al., 1996;Stolle et al., 2009).In our analyses, the universal form (determined from Eq. ( 2) based on the estimated α and C 1 ) fits the empirical K(q) (determined from the logarithmic slopes of trace moments) quite well.A multifractal phase transition is observable in the plot (Fig. 4) of empirical and theoretical K(q) curves, indicating that the measured moments will only have the theoretical K(q) for q below a critical moment q c .Beyond q c there is a multifractal phase transition where K(q) becomes asymptotically linear.The linear behavior of the empirical moment scaling function is either due to sampling limitations (i.e., second-order multifractal phase transition; Schertzer and Lovejoy, 1992) or its association with a divergence of statistical moments (i.e., firstorder multifractal phase transition; Schertzer and Lovejoy, 1992).In the first-order multifractal phase transition, q c corresponds specifically to maximum singularity measured and is associated with the occurrence of very rare and violent singularities, whereas in the case of a second-order multifractal phase transition, q c corresponds to the maximum singularity effectively measurable from a finite sampling (Seuront et al., 1999).
In order to differentiate between first-and second-order multifractal phase transitions, we compare the theoretical value of the critical moment q s with the empirical critical moment q c calculated from the K(q) curve.The theoretical value of the critical moment q s can be computed as (Schertzer and Lovejoy, 1992): If q c ≈ q s , the phase transition can be termed as a secondorder multifractal phase transition wherein the critical moments are only related to sampling limitations.Also, if q c < q s , the critical moments q c are independent of the sampling and characterize the occurrence of very rare and violent singularities present in the data set (i.e., first-order multifractal phase transition) (Seuront et al., 1999).Using the values of C 1 and α, the average q s values computed for coarse sediment region are found to be 2.32 and 2.16, respectively, at 33 and 210 kHz, whereas in fine sediment region, the average q s values are observed to be 4.22 and 3.97, respectively, at 33 and 210 kHz.The estimated q s values indicate a first-order multifractal phase transition (q c < q s ) in the dual-frequency echo envelopes with the occurrence of rare and violent singularities in the data set.The detection of the presence of a first-order multifractal phase transition possibly suggests that the time-dependent dual-frequency seafloor backscattering could be a self-organized critical (SOC) process.q -2.0 1.2 q -0.5 q -2.0 q -1.2 q -0.5 q -2.0 q -1.2 q -0.5 q -2.0 q -1.2 q -0.5 L = 5 ms L = 2.8 ms The q values of the each trace moments are varied between 0 and 2, in intervals of 0.1 (from top to bottom, q = 2-1, 0.1, 0.9, 0.2, 0.8, 0.3, 0.7, 0.4, 0.6 and 0.5).The multiscaling of trace moments generally holds quite well up to q < 1.6.The deviation of trace moments (dots) from the fitted linear curve is indicative of "break" in the multiscaling.The break is less apparent for low q(< 1.6) values, but becomes conspicuous for large q(> 1.6) values.The break associated with the time-dependent seafloor backscattering is possibly due to the collective effect of inherent heterogeneities present in the seafloor (mainly because of the coarse sand particles, shells, gas bubbles, benthic organisms and sediment layers).Depending on acoustic wavelength/frequency, the individual features such as shells and other roughness elements at the sediment-water interface may be more appropriately characterized as discrete scatters than as micro topography (Jackson and Richardson, 2007).An appropriate assessment in this regard is difficult due to lack of supporting data (similar results were also reported earlier by Gagnon et al. (2006) and Lovejoy and Schertzer (2007b) while analyzing the high-resolution Lower Saxony DEM data over Germany.The break observed in their analyses was due to the effect of "trees" evident in the high-resolution topography data.).

The universal multifractal parameters
The dual-frequency universal multifractal parameters were analyzed along with the ground truth values of the mean grain size (M ϕ ) to understand their relationships.The C 1 values are negatively correlated with the measured M ϕ , having correlation coefficients of −0.96 and −0.92, respectively, at 33 and 210 kHz (Fig. 5).The C 1 values decrease with increasing weight percentage of both silt and clay fraction (with M ϕ > 4).Concurrently, the C 1 values increase linearly with increasing percentage content of sand fraction (with M ϕ < 4).The range of C 1 values associated with coarse and fine sediments vary between 0.171-0.249and 0.035-0.089and between 0.180-0.294and 0.051-0.090,respectively, at 33 and 210 kHz.Moreover, in the coarse sediment region, the average C 1 values are restricted to values around 0.203 ± 0.0222 and 0.209 ± 0.0372, respectively, at 33 and 210 kHz.In the fine sediment region, the average C 1 values are found to be within 0.066 ± 0.0162 and 0.0722 ±v0.0115, respectively, at 33 and 210 kHz.In brief, the C 1 values are well clustered at both the acoustic frequencies, with fewer fluctuations for the fine sediment as compared with the coarse sediment region found at deeper depths (Fig. 5).The relatively low C 1 values attributed to the fine sediment region indicate that the field values (pressure values) are close to the mean values.The α values of the dual-frequency echo envelopes show identical trends (≈ 1.93), expressing similar degrees of multifractality in coarse and fine sediment provinces (Figs. 6 and  7).The H values signifying the degree of smoothness of the data associated with coarse and fine sediments vary between 0.857-0.999and 0.888-0.984and between 0.887-0.986and 0.804-0.913,respectively, at 33 and 210 kHz.The scatter diagram (Fig. 7) of H values at 33 and 210 kHz reveals Comparison between the empirical (solid curve) and theoretical (dash curve) moment scaling function K(q) as a function of q.The empirical K(q) values of the representative clayey-silt and sand substrates are determined from the logarithmic slopes of trace moments illustrated in Fig. 3.The theoretical K(q) values are calculated using Eq. ( 2) based on the computed α and C 1 (shown in the figure).The linear deviation of theoretical curves beyond q ≥ q c > 1.25 is indicative of a first-order multifractal phase transition, caused by the rare and violent singularities in the dual-frequency data.
that the values computed at 210 kHz are confined within 0.874 ± 0.0406 in fine sediments and within 0.951 ± 0.0268 in coarse sediments.However, at 33 kHz the calculated H values do not exhibit any apparent trend distinguishing between the fine and coarse sediment provinces.The obtained results provide a construal that is similar to the previously estimated sediment geoacoustic inversion results (discussed in the Sect.3.4).Generally, the dual-frequency universal multifractal parameters among the coarse and fine sediments show subtle differences in α and H , whereas the codimension parameter C 1 representing the sparseness of the data varies (Fig. 7 and Table 1).This suggests that the physics of the scattering mechanism (discussed in Sect.3.3.1)responsible for the variation in C 1 is different.

Possible influence of the seafloor backscattering process
The backscattering from the seabed can be generally attributed to two contributing factors, namely interface and volume scattering (Sternlicht and de Moustier, 2003a).The strength of the backscattered signal is primarily controlled by the acoustic frequency, the acoustic impedance contrast between water and sediment, the contributions from seafloor interface roughness, as well as the sediment volume heterogeneity.The interface scattering is governed by the microscale roughness of the seafloor facets coupled with the acoustic impedance.A part of the transmitted acoustic energy penetrates the sediments and is reflected back by the volume heterogeneities.Such scattering mechanism is normally referred to as volume scattering (mainly due to the coarse sand particles, shells, gas bubbles, benthic organisms and sediment layers).
The shape of the echo envelope has two distinct parts, the initial part and the tail part.The initial part of the data represents the reflection from the water-sediment interface (interface scattering), and the tail portion corresponds to the backscatter from the sediment volume (volume scattering).A significant contribution, due to the volume scattering (in addition to the interface scattering) from various scatterers, is expected to be dominant for acoustically soft sediments such as mixtures of clayey silt and silt.Accordingly, the contribution of sub-bottom scattering becomes prominent near the tail portion of the echo envelope from the soft sediments (De and Chakraborty, 2011).Besides, acoustically soft sediments are penetrated more deeply by the acoustic signal and produce longer and corrugated echoes than hard sediments.The scattering process takes place exclusively at the surface for  acoustically hard sediments, i.e., mixtures of silty-sand and sand sediments.The scattering processes described herein determine the statistical and geometrical properties of the data, resulting in the variation in the estimated C 1 parameter.

Relationship with fractal dimensions
Several studies while comparing the fractal dimension of the echo envelopes with the ground truth sediment have concluded that the fractal dimension (as a measure of complexity and roughness) is a good descriptor of a bottom type in the investigated area (Tegowski and Lubniewski, 2000;Tegowski et al., 2003;Tegowski, 2005;van Walree et al., 2005;De and Chakraborty, 2009).The fractal dimension describes the statistical and geometrical properties of the data.Variations in the fractal dimension of the echo envelope carry information concerning the fractal structure of the sediment layers signifying the hardness of the seabed.The fractal dimension of the echo envelopes reflected from the fine seafloor has been found to be higher as compared with the coarse sediment region (Tegowski and Lubniewski, 2000).The low values of fractal dimension in the coarse region have been controlled by the dominant interface scattering due to limited bottom penetration of the acoustic signals.Besides, the fine region reflects the environment with deeper acoustic penetration and the fractal structure of the sediment layers are more apparent in the recorded echo, indicating higher fractal dimension values.
The C 1 parameter represents the codimension of the set of points that gives the dominant contribution to the mean of the conserved multifractal flux.The corresponding fractal dimension can be expressed as d − C 1 , where d denotes the standard dimension of the space (Lovejoy et al., 2001).The relatively high C 1 values correspond to a very sparse process.The field values contributing to the mean behavior in such an instance are violent and confined to a very sparse set signifying low fractal dimension.Conversely, a low C 1 implies a more uniform and less extreme process.Appropriately, the relatively low values of C 1 attributed to the fine sediment region (with higher fractal dimension) indicate that the field values (pressure values) are close to the mean values as compared with the coarse sediment region (Fig. 5).
The stochastic-based multifractal analysis followed herein has several advantages over standard statistical approaches,  as it characterizes the local scale properties of the data in addition to its global properties.Correspondingly, it is possible to quantify the statistical distribution of the local singularities (i.e., local multifractal exponents) present in the data (see Lovejoy et al., 2009).

Comparison with inversion results
The acoustic backscatter data obtained from the echosounding systems can be matched with theoretical scattering models to interpret the information embedded in the data (Jackson et al., 1986).The numerical approach employed for extracting information from the data is commonly referred to as "inversion modeling".The inversion modeling primarily involves physics-based inversion of echosounding data to obtain the upper-layer seafloor roughness parameters, namely the sediment mean grain size (M ϕ ); spectral parameters at the water-seafloor interface (γ 2 , w 2 ); and sediment volume parameter (σ 2 ), which can be used to examine the fine scale seafloor processes (Sternlicht and de Moustier, 2003a, b;De and Chakraborty, 2011;Haris et al., 2011).
The seafloor "roughness power spectrum" estimated from the echo data characterizes the size and periodicity of the seafloor height fluctuations as a function of the spatial fre-quency.The roughness power spectrum is often parameterized using a power law by slope and intercept of a linear regression line through the points of the periodogram estimate in log-log space.Indeed, the parameters γ 2 and w 2 used in the scattering models are the slope and intercept, respectively, of the 2-D roughness power spectrum, which are estimated from the 1-D power-law values (Haris et al., 2011).A wide range of 2-D roughness power spectrum parameters of the study area are available (De and Chakraborty, 2011) and offer an opportunity to determine their relationship with the presently estimated universal multifractal parameters.
Parameter H describing the height statistics of the data has been empirically determined from the spectral slope β and K(2) using the relation β =1+2H -K(2) (Table 1).The differences in H principally reflect variations in the spectral slopes (Lovejoy et al., 2001) (and thereby the corresponding γ 2 values).With reference to the inversion modeling study carried out by De and Chakraborty (2011), in the coarse sediment region, the average γ 2 values were restricted to values around 3.23 ± 0.071 and 3.16 ± 0.047, respectively, at 33 and 210 kHz.In the fine sediment region, the average γ 2 values were found to be within 3.22 ± 0.074 and 3.30 ± 0.037, respectively, for 33 and 210 kHz.The scatter diagram (Fig. 3b    of De and Chakraborty, 2011) between the estimated mean values of M ϕ and γ 2 at 210 kHz indicated that the values of γ 2 were confined within 3.21-3.4 in fine sediments and within 3.0-3.21 in coarse sediments.In contrast, the estimated mean values of γ 2 at 33 kHz inversions did not exhibit any apparent trend to discriminate between the fine and coarse sediment provinces.The subtle difference in the computed γ 2 (or the spectral slope β) between coarse and fine sediments con-form to the meager variation in the computed H parameter (Fig. 7), particularly at 33 kHz.The scatter diagram between the estimated mean values of M ϕ and w 2 (Fig. 3a of De and Chakraborty, 2011) revealed that, in the coarse sediment region, the average w 2 values were restricted to values around 0.00356 ± 0.00047 and 0.00365 ± 0.00101, respectively, at 33 and 210 kHz, but that in fine sediments, the average w 2 values were found to vary between 0.000461 ± 0.00013 and 0.000605 ± 0.000042, respectively, at 33 and 210 kHz.The computed w 2 (or intercept) values were well clustered at both the acoustic frequencies, having fewer fluctuations for the fine sediment as compared with the coarse sediment region.Likewise, the parameters C 1 revealing the noise statistics of the data are well clustered at both the acoustic frequencies with fewer fluctuations for the fine sediment as compared with the coarse sediment region.It is observed that the relatively higher values of w 2 and C 1 are associated with coarse sediments, while the lower values of w 2 and C 1 are the characteristics of fine sediments (Fig. 8).
As pointed out in the introduction, the model-based methods can help interpret the echo signal of the seafloor sediment properties (M ϕ ) and micro roughness parameters (γ 2 , w 2 ).However, the calculation of a correct set of geoacoustic parameters gets convoluted by the large number of good fits existing in the multidimensional search space.Accordingly, it is possible to obtain convincing model-data fits in the search space that do not necessarily correspond to the correct set of geoacoustic parameters (Sternlicht and de Moustier, 2003b).Moreover, the physics-based models are valid only for a certain range of frequencies and sediment types (Amiri-Simkooei et al., 2011), such that the direct inversion of an acoustic signal is unlikely without setting the limits of geoacoustic parameters for a known seabed sediment.In contrast, the statistically based empirical methods rely on the analyses of certain echo signal features that are correlated with sediment properties.These methods are relatively easy to implement in view of the computation time involved.However, proper ground-truth measurements are imperative to validate and interpret the results.

Concluding remarks
The scaling and multifractality of the SBES echo envelopes (at 33 and 210 kHz) have been demonstrated using a stochastic-based universal multifractal framework.The variations in the three fundamental scale invariant parameters (α, C 1 , and H ) and their behavior with respect to the ground truth sediment information can delineate different bottom types in the investigated area.The evidence for the first-order multifractal phase transition (with the divergence of higher order statistical moments) further reveals the hierarchy of multiplicative cascade dynamics associated with the echo envelopes.
The universal multifractal parameters among the coarse and fine sediments show subtle differences in α and H , whereas the codimension parameter C 1 representing the sparseness of the data varies.The C 1 values are well clustered at both the acoustic frequencies, demarcating the coarse and fine sediment provinces.The computed α values show an identical trend (≈ 1.93), expressing a similar degree of multifractality in the coarse and fine sediment regions.The minute variation in the H parameter ascribed to the coarse and fine sediment locations is well corroborated by the previously estimated sediment geoacoustic inversion results.In the context of multifractal analyses the 210 kHz appears to be marginally better as compared with 33 kHz.This could be due to the fact that the lower acoustic frequencies penetrate relatively more into the substrates, whereas higher frequencies have a better resolving capability.The final outcome of the multifractal analyses underpins the hitherto applied model-based seafloor characterization and helps foster research on the empirical method-based feature discrimination.
Fig. 1.(a) GIS-based sediment distribution map of the study area showing the acoustic and sediment data acquisition locations.The acoustic data were sampled with RESON-NS 420 dual-frequency SBES along the three tracks.The ground truth sediment information was collected using a Van Veen grab.Panel (b) represents the graphical abstract of the methodology implemented in the study.The dual-frequency echo envelopes (red curve) from 20 locations are subject to stochastic-based universal multifractal analyses to verify the existence of multiscaling.The shaded region represents the variability of the data.

Fig. 2 .
Fig. 2. Mesh plot illustrating the effect of echo alignment technique.The initial rise time and the time of peak amplitude of the raw backscatter envelope records are selected as the aligning features.The resulting aligned echo envelopes converted to the pressure values were ensemble averaged to obtain a representative stable acoustic signal (at each location) prior to multifractal analyses.

Fig. 3 .
Fig. 3. Panels (a) and (b) shows the log/log plot of the normalized trace moments as a function of the scale ratio λ = L/ l at 33 and 210 kHz.The multiscaling of the statistical moments of the representative clayey silt and sand substrate is verified by the straightness of the curves.The q values of the each trace moments are varied between 0 and 2, in intervals of 0.1 (from top to bottom, q = 2-1, 0.1, 0.9, 0.2, 0.8, 0.3, 0.7, 0.4, 0.6 and 0.5).The multiscaling of trace moments generally holds quite well up to q < 1.6.The deviation of trace moments (dots) from the fitted linear curve is indicative of "break" in the multiscaling.The break is less apparent for low q(< 1.6) values, but becomes conspicuous for large q(> 1.6) values.The break associated with the time-dependent seafloor backscattering is possibly due to the collective effect of inherent heterogeneities present in the seafloor (mainly because of the coarse sand particles, shells, gas bubbles, benthic organisms and sediment layers).Depending on acoustic wavelength/frequency, the individual features such as shells and other roughness elements at the sediment-water interface may be more appropriately characterized as discrete scatters than as micro topography(Jackson and Richardson, 2007).An appropriate assessment in this regard is difficult due to lack of supporting data (similar results were also reported earlier byGagnon et al. (2006) andLovejoy and Schertzer (2007b) while analyzing the high-resolution Lower Saxony DEM data over Germany.The break observed in their analyses was due to the effect of "trees" evident in the high-resolution topography data.).
Fig. 4.Comparison between the empirical (solid curve) and theoretical (dash curve) moment scaling function K(q) as a function of q.The empirical K(q) values of the representative clayey-silt and sand substrates are determined from the logarithmic slopes of trace moments illustrated in Fig.3.The theoretical K(q) values are calculated using Eq.(2) based on the computed α and C 1 (shown in the figure).The linear deviation of theoretical curves beyond q ≥ q c > 1.25 is indicative of a first-order multifractal phase transition, caused by the rare and violent singularities in the dual-frequency data.

Fig. 5 .
Fig. 5. Frequency-wise scatter plot of the measured M ϕ and computed codimension parameter C 1 at 33 and 210 kHz.The shaded region demarcates the boundaries of C 1 in coarse and fine sediment provinces.

Fig. 6 .
Fig. 6.A quantitative comparison between the scale invariant multifractal parameters (α, C 1 , and H ) at 33 and 210 kHz.The shaded region represents the coarse and fine sediment provinces.

Fig. 7 .
Fig. 7. Panel (a) shows the comparison of the computed multifractal parameters (α, C 1 , and H ) at frequencies 33 and 210 kHz.The C 1 values are well clustered, demarcating the coarse and fine sediment provinces.The horizontal dash line (red) demarcates the H values at 210 kHz in fine and coarse sediment provinces.The variations in multifractal parameters for different frequencies with the same sediment types at each location are represented in panel (b), the shaded region corresponds to coarse sediment locations.

Fig. 8 .
Fig. 8. 3-D scatter diagram showing clustering among M ϕ , C 1 and w 2 at 33 and 210 kHz.Regions comprising coarse and fine sediments can be seen clearly delineated.

Table 1 .
Summary of universal multifractal parameters.