A Method for Estimating Dominant Acoustic Backscatter Mechanism of Water-Seabed Interface via Relative Entropy Estimation

It is important to distinguish the dominant mechanism of seabed acoustic scattering for the quantitative inversion of seabed parameters. An identification scheme is proposed based on Bayesian inversion with the relative entropy used to estimate dominant acoustic backscatter mechanism. DiffeRential Evolution Adaptive Metropolis is used to obtain samples from posterior probability density in Bayesian inversion. Three mechanisms for seabed scattering are considered: scattering from a rough water-seabed interface, scattering from volume heterogeneities, andmixed scattering from both interface roughness and volume heterogeneities. Roughness scattering and volume scattering are modelled based on Fluid Theories using Small-Slope Approximation and SmallPerturbation Fluid Approximation, respectively. The identification scheme is applied to three simulated observation data sets. The results indicate that the scheme is promising and appears capable of distinguishing sediment volume from interface roughness scattering and can correctly identify the dominant acoustic backscatter mechanism.


Introduction
Acoustic scattering of water-seabed interface is the main source of underwater reverberation, which could affect the performance of sonar system when it is used to detect and classify underwater targets [1].If we are interested in the seabed, the acoustic scattering is the main way to obtain information of seabed properties.It is believed that seabed scattering consists mainly of two parts: interface rough scattering and inhomogeneous sediment volume scattering [2].In order to mitigate the reduction of sonar performance and enhance the acquisition of seabed properties, it is necessary to identify the dominant acoustic scattering mechanism and choose an accurate model of seabed scattering at the location of interest [3].This paper develops a Bayesian approach to identify dominant acoustic backscatter mechanism via relative entropy estimation, figuring out whether interface scattering, volume scattering, or both interface and volume scattering is the main factor.
Scattering mechanisms are distinguished using the relative entropy [4,5], also known as Kullback-Leibler divergence, which is an optimal criterion used to describe the difference between two probability distributions (prior probability distribution and posterior probability density).Compared to other optimal criteria, the relative entropy is suitable for nonlinear model and non-Gaussian distribution of parameters, measuring the amount of information generated by observations and selected model.The motivation of this paper is the application that relative entropy is used to measure information provided by an experiment design [6].In this paper, the design of the numerical experiment, such as grazing angle and acoustic frequency of the scattering strength measurement, is fixed while different models are used.Then different relative entropies correspond to different models and different models correspond to different main scattering mechanisms.In order to obtain an accurate estimation of relative entropy, the main task is to solve the posterior probability density (PPD) which can be calculated by Model-based Bayesian geoacoustic inversion.
Since the start of geoacoustic inversion more than two decades ago, there are many geoacoustic bottom interaction models, including roughness scattering models and sediment volume scattering models based on different media theories (such as Fluid Theories, Elastic Theories, and Poroelastic Theories), that are established and tested [7,8].Although different models differ in the fit degree with the measured data, these models prove their effectiveness with varying degrees.Specifically, considering that the purpose of this paper is to distinguish the different scattering mechanisms, a simplified forward model which takes sediment as fluid is used to predict scattering at the water-sediment interface and/or within the volume of sediment.The seafloor roughness and sediment heterogeneities are taken as random processes which follow the pure power law.Small-Slope Approximation and Small-Perturbation Fluid Approximation are used to calculate the interface scattering and volume scattering, respectively.
Bayesian inference has been widely used in geoacoustic inversion to obtain parameter uncertainty analysis instead of point estimation based on global optimization techniques.The key task in Bayesian inference is to solve the posterior distribution.However, geoacoustic bottom interaction models are strongly nonlinear that means the posterior distribution has no analytical form.Monte Carlo simulation methods can be used to generate a sample from the posterior distribution, also referred to as the target distribution.Then the samples can be used to approximate the posterior distribution.Many iterative methods are used to generate samples that are as close to the target distribution as possible.Most of these methods are based on Markov Chain Monte Carlo (MCMC) approach, which exploit a Markov chain that generates a random walk through the parameter search space to obtain solutions with stable frequencies stemming from a stationary distribution.The earliest MCMC approach is the random walk Metropolis.Considering nonsymmetrical jumping distributions, Hastings extended the random walk Metropolis algorithm to more general cases, which is well known as the Metropolis-Hastings (M-H) algorithm [9].In order to increase sampling efficiency of complex posterior distributions involving long tails, correlated parameters, multimodality, numerous local optima, more and more improved algorithms have been proposed such as adaptive proposal (AP) [10], adaptive Metropolis (AM) [11], delayed rejection adaptive Metropolis (DRAM) [12], differential evolution Markov chain (DE-MC) [13], Fast Gibbs sampling (FGS) [14,15], and so on.The sampling method used in this paper is a relatively efficient one entitled DiffeRential Evolution Adaptive Metropolis (DREAM) [16], which can maintain detailed balance and ergodicity with excellent performance on a wide range of problems involving nonlinearity, high dimensionality, and multimodality.
In this paper, the high-frequency seafloor acoustics models [17] are adopted, and it is assumed that the sound wave only interacts with the surface of seafloor without considering the influence of the stratification of seabed sediments.Three simulated data sets (two representing a muddy sediment and one representing a sandy sediment) are used in the Although the selected seabed parameters and simulated data sets may differ from the actual seabed situation, it is enough to verify the effectiveness of proposed method.The remainder of this paper is organized as follows.Section 2 describes the scatter models for calculation of scattering strength, Bayesian inference, and sampling method.Section 3 presents the simulation study of the method to identify the dominant backscatter mechanism.Finally, Section 4 summarizes and discusses this work.

Materials and Methods
For each set of observation data, the relative entropy estimation steps of identification scheme are as follows: (a) parameter inversions are carried out based on the assumption of different dominant backscatter mechanisms; (b) the sample values from the posterior probability density are used to calculate the relative entropy; (c) the values of relative entropy based on different assumptions are compared with each other to identify the dominant backscatter mechanism.To implement the above steps, this section describes the calculation process of scatter strength, Bayesian inference, and sampling method.

Scatter Model.
Although the separation of scattering into separate, independent roughness and volume components is somewhat artificial [17], it is meaningful to understand the scattering mechanism better by discussing them separately.Before introducing the model calculation, the geometric model of sound propagation is established as shown in Figure 1.
The in-water wave vectors of incident and scattered directions are defined as follows: where   = 2/  is the acoustic wavenumber in water.  and  are velocity of sound in water and acoustic frequency.
Thus, the horizontal components and vertical component of the above wave vectors, respectively, are where the uppercase, boldface symbols represent the horizontal components.The following wave vector differences needed in the model calculation are defined: Based on Snell's law, the wave vectors for the incident and scattered compressional waves in the sediment are where   is the acoustic complex speed ratio and its expression will be given in the next part.Then the in-sediment three-dimensional Bragg vector is In the calculation below, the magnitude of a difference vector is denoted without boldface.Assuming that seafloor is isotropic in the statistical sense, we may set 0  = 0  = 0,   =  −   for backscatter.

Roughness Scattering Model.
There are two most widely used approximations for acoustic scattering by seafloor roughness, which are the small-roughness perturbation method [18] (sometimes known as Rayleigh-Rice perturbation theory) and the Kirchhoff approximation [19] (also known as the tangent-plane approximation).The perturbation theory tends to be most accurate for scattering at wide angles relative to the specular (flat-interface reflection) direction and the Kirchhoff approximation is better for scattering near the specular direction.The complementary behaviour of the Kirchhoff and perturbation approximations has led numerous investigators to combine the two in the composite-roughness approximation.Luckily, there is a kind of approximation method entitled small slope approximation [20,21] which can provide a single expression that covers all angles and is likely to be at least as accurate as either the Kirchhoff or perturbation approximations.To predict roughness scattering strength, five parameters are needed, which are roughness spectral exponent  2 , roughness spectral strength  2 , ratio V  of sound velocity in sediment to   , ratio   of sediment density to water density and, ratio   of the imaginary and real parts of complex compressional speed.The final expression of scattering strength   and calculation process are as follows [17]: Assuming that the roughness follows pure power law, the Kirchhoff integral   is where J 0 is the zeroth-order cylindrical Bessel function of the first kind and  and  are as follows: 2 ℎ is the square of the "structure constant" which is related to the parameters of the power-law spectrum through where Γ is the gamma function.The final expression of factor   is where   is the flat-interface reflection coefficient, where   is the acoustic complex speed ratio which is defined as So far, the scattering strength   can be finally obtained.

Volume Scattering Model.
Compared to the hard sediment, the soft sediment with heterogeneities could produce appreciable levels of scattering strength without considering the possible presence of embedded shells and bubbles.The reason is that the sound speed in soft sediment is similar with   and the sound can easily transmit into and out of sediment.Typical soft sediment is mud compared with sand.Six parameters are needed for the calculation of the volume scattering strength.V  ,   , and   have been defined before.The other three parameters are density fluctuations spectral exponent  3 , density fluctuations spectral strength  3 , and  which is a dimensionless parameter describing the correlation between the spectrums of sediment compressibility fluctuations and density fluctuations.If the volume scattering strength is independent of depth, the expression for  V is [22,23] ) . ( Equation ( 21) forms the basis for several sediment volume scattering models, which differ only in the assumption and approximation used to obtain the volume scattering cross section  V .Here the Small-Perturbation Fluid Approximation, which is also called Born approximation, is used to calculate where   is the spectrum for density fluctuations with the pure power law form: At last, the scattering strength of roughness and volume scattering (mixed scattering model) can be obtained through ) . (24)

Bayesian Inference for Relative Entropy.
A more detailed description of Bayesian inference applied to different inverse problems can be referred to [25][26][27].Bayesian inference is an ideal tool to handle uncertainties analysis instead of point estimation based on global optimization algorithm.Intuitively, all model parameters are considered as random variables.When a set of observations (scattering strength)  and the scattering model are determined, the PPD of parameter vector  can be written as where () is the prior probability distribution of  and ( | ) is conditional probability distribution of observations and is generally considered as likelihood function.() is independent with  and is usually considered as constant.
Then the PPD can be written as Assuming zero-mean Gaussian-distributed errors, the likelihood can be expressed as where   is the dimension of , () denotes the vector predicted by the scattering model, and   is a diagonal data covariance matrix.One-and two-dimensional marginal PPD can be obtained as follows: where the non-boldface symbols   ,   represent members of vector .For the prior probability distribution and PPD, the expression of relative entropy is This quantity is nonnegative and nonsymmetric and reflects the difference of information carried by the two distributions.The larger relative entropy means that the posterior probability is more concentrated compared to the prior probability.For the inversion results in this paper, the more concentrated posterior probability means the observation data and the model match better and corresponds to the dominant backscatter mechanism.The integral equations ( 28), (29), and (30) that need to be solved are rewritten in the following uniform form Based on the theory of Monte Carlo, if the integral does not have an analytic solution, we can rewrite (31) as follows: So the integral becomes the mean of   () on the distribution function ().If we get a large number of samples  1 , 2 . ..  from (), the integral can be approximately equal to If samples are generated from ( | ), the calculation process can be further simplified by In some literatures, multiple maximum posterior (MAP) estimation methods [28,29] are used to get samples with unknown () such as genetic algorithm (GA) and obtain parameter uncertainty estimation by equation (34).The samples may be too concentrated near the MAP estimation and do not represent the true target distribution.So we need more rigorous sampling method which will be described in the next part.

DiffeRential Evolution Adaptive Metropolis.
In order to obtain samples from target distribution, MCMC algorithm is the main and effective method.Its basic idea is constructing a Markov chain which can explore the parameter space and gradually converge to its equilibrium distribution.The key that affects the efficiency of the algorithm is whether an optimal starting point and proposal distribution can be found.The starting point can be a MAP estimation which we get from global optimization algorithm such as GA, Simulated Annealing (SA).Compared with starting point, proposed distribution has a larger impact on convergence efficiency.The proposed distribution actually affects the parameter perturbation and the efficiency of exploring the entire parameter space.Compared to many other improved MCMC algorithms, DREAM has proved its excellent performance on target distribution sampling involving nonlinearity, high dimensionality, and multimodality [16].
In general, the probability to accept candidate state   as the next state of  −1 within Markov chain is where

Results and Discussion
In this paper, the dominant acoustic backscatter mechanism is identified by the relative entropy.For each set of observations, three inversions are conducted based on the assumption of different dominant acoustic backscatter mechanisms, which are interface roughness scattering, volume scattering, and both interface roughness and volume scattering (mixed scattering).Once all inversions are completed, the relative entropy for each inversion is calculated to measure the amount of information generated by observations and selected model.As described in Section 2.2, when the observation data and the used model match best, the largest relative entropy which corresponds to the main scattering mechanism can be obtained.The modification process of acoustic scattering model motivated by the mismatch of data and model is actually a process of modification and addition of additional physical mechanisms [8].Models that consider more scattering mechanisms will match the data best.In other words, the inversion results based on the assumption of mixed scattering may have the largest relative entropy.That doesn't mean the dominant acoustic backscatter mechanism is mixed scattering due to the good fitting performance of mixed scattering model.We also have to consider whose relative entropy of the assumption is very close to that of mixed scattering.Only if the relative entropys of the other two mechanisms are both smaller than that of the mixed scattering, it can be considered that the dominant backscatter mechanism is mixed scattering.In this section, the proposed identification method is validated by nine simulated-data inversions.There are three sets of geoacoustic parameters, which are selected referring to [17], representing a sandy seabed and two types of muddy seabed.Parameter set 1 represents the sandy seabed corresponding to the dominant backscatter mechanism of roughness scattering.Parameter set 2 represents the muddy seabed with dominant backscatter mechanism of volume scattering.The last parameter set of muddy seabed corresponds to the dominant backscatter mechanism of mixed scattering.True values of the parameter sets needed for the calculation of scattering strength are given in Table 1.For each sediment type, the scattering strength is computed for roughness scattering, volume scattering, and mixed scattering.Random errors are added to the true predicted scattering strength through where () denotes the vector of true predicted scattering strength and  denotes the observation data vector, which is used for inversions.Each element of  is drawn from the standard normal distribution."∘" represents the operation: ( ∘ ) , = () , () , .The angular and frequency ranges for () are { 10 ∘ ∼85 ∘ } and {20 kHz, 40 kHz, 60 kHz}, respectively.Each of the three observation data sets is inverted once assuming each of the three types of scattering.
Figure 2 shows the marginal PPD of each parameter from the inversion results for the true roughness scattering In addition to the uncertainty of the parameters, we can also obtain the uncertainty of predicted scattering strength calculated by the uncertainty of the parameters, as shown in Figure 3.It can be seen that the trend of uncertainty is consistent with the true scattering strength.Compared with large (  > 60 ∘ ) and small (  < 20 ∘ ) grazing angle, the uncertainty of predicted scattering strength for intermediate range of grazing angle is larger due to the larger error of observations.The drop above the critical angle of 30 ∘ is due to the increased penetration into sediment and it is a unique characteristic of sandy seabed [30].
The inversion results in Figure 4 for true volume scattering are similar with the results in Figure 2. difference is that inversion results of volume scattering model are better than that of roughness scattering model.Due to the real dominant backscatter mechanism, parameters  2 ,  2 , which describe the interface roughness, are not properly estimated.The inversion result of   of mixed scattering model is obviously better than that of volume scattering model and roughness scattering model.For the other five parameters, the PPDs of mixed scattering model and volume scattering model are similar in general.
There is no drop for the predicted scattering strength of muddy seabed in Figure 5 since mud has a similar sound speed with water.This feature is often used to distinguish sandy and muddy seabed.Volume scattering strength does not keep increasing at nearly vertical angles.In a mixed scattering situation, the main contribution of volume scattering to total scattering focuses on the small to medium grazing angle and sometimes eliminates the drop above the critical angle of rough scattering as shown in Figure 7.
Figure 6 shows the inversion results for the observations of scattering strength with true mixed scattering.Obviously, the inversion results of mixed scattering model are more satisfied than that of roughness and volume scattering model.However, it can be seen that the PPDs of parameters  3 and  3 are relatively poor compared with the results in Figure 4.
Each of the three dada sets is inverted once assuming each of the three types of scattering, and the relative entropy is calculated for each of the nine inversions which are shown in Table 2.It is not difficult for us to find that the magnitude of the relative entropy corresponds essentially to the quality of the previous inversion results.For the true roughness scattering, the relative entropy of assumed roughness scattering is very close to that of assumed mixed scattering which is the maximum.The other two largest relative entropys' assumptions all correspond to the true scattering mechanisms.The relative entropy correctly identifies the dominant backscatter mechanism in all cases.For the true mixed scattering, the relative entropy of assumed volume scattering is far less than that of assumed roughness scattering, which means the contribution of roughness scattering to total scattering is much larger.This difference also explains why the PPDs of parameters  3 and  3 in Figure 6 are relatively poor compared with the results in Figure 4.

Conclusions
This paper proposed a method to remotely identify the dominant acoustic backscatter mechanism of water-seabed interface based on the relevant parameters uncertainty of seabed sediments, respectively.An identification scheme based on the relative entropy is proposed for determining the dominant acoustic backscatter mechanism.The relative entropy is an optimal criterion used to describe the difference between prior probability distribution and PPD and it is suitable for nonlinear model, measuring the amount of information generated by observations and selected model.The identification scheme is tested using inversions of three simulated scattering data sets and it correctly identifies the dominant acoustic backscatter mechanisms.The simulation results indicate that the proposed method is promising and appears capable of distinguishing sediment volume from interface roughness scattering.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Figure 2 :
Figure 2: Marginal PPD of each parameter assuming each of the three types of scattering (Parameter set 1).

Figure 3 :
Figure 3: True predicted scattering strength ( * ), observations of scattering strength (∘) and uncertainty (color part) of predicted scattering strength by inversion results of mixed scattering model (Parameter set 1).

Figure 4 :Figure 5 :
Figure 4: Marginal PPD of each parameter assuming each of the three types of scattering (Parameter set 2).

Figure 6 :Figure 7 :
Figure 6: Marginal PPD of each parameter assuming each of the three types of scattering (Parameter set 3).
The candidate state   will be accepted as the next state of  −1 if  ≤   .Otherwise, the next state   will remain the same with  −1 .The efficiency of this seemingly simple process will be significantly affected by parameter perturbation until the chain finally explores the entire parameter space.DREAM uses N Markov chains and differential evolution to expedite the exploration process.The candidate state of the No. ( = 1, 2 . . ., ) chain is updated by    =    + (1   +    )  (,  )   parameters changed by    and keep the rest   −   parameters unchanged within the No.i chain.The values of vector    and    are sampled independently from a multivariate normal distribution    (0, ) and a multivariate uniform distribution    (−  ,   ), respectively, with default values   = 0.1 and  = 10 −12 . denotes the number of chains used to generate the perturbation.  and   denote the No.j element of vectors  and  consisting of  integers drawn without replacement from {1, 2, . . ., }. (,  ) = 2.38/√2   is the jump rate which is periodically with a 20% chance to enhance the probability of jumps between disconnected modes of the target distribution.Then the candidate state of No.  chain at  is (  ) and ( −1 ) are the posterior probability density of parameter vectors at different time and (  →  −1 ) and ( −1 →   ) are the conditional probabilities of chain moves from   to  −1 and  −1 to   , respectively.(−1→   ) = min [1,  (  )  ( −1 )] .(36)Then a random sample  taken from a uniform distribution (0, 1) will be compared with   ( −1 →   ).

Table 1 :
True parameter values for simulated inversions.

Table 2 :
Relative entropy for different datasets based on different assumptions.Row names are the assumed backscatter mechanisms, column names correspond to the true backscatter mechanisms.Bold values indicate largest relative entropy in the group.