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

Authors: Bo Zou; Jingsheng Zhai; Jian Xu; Zhaoxing Li; Sunpei Gao
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

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, and mixed scattering from both interface roughness and volume heterogeneities. Roughness scattering and volume scattering are modelled based on Fluid Theories using Small-Slope Approximation and Small-Perturbation 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.


---

## Body

## 1. 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 study of dominant acoustic backscatter mechanism estimation. The dominant scatter mechanisms for the three data sets are volume scattering, interface scattering, and mixed scattering, respectively. Based on the assumption of different dominant backscatter mechanisms, inversions are carried out for each data set. The estimation results can correctly identify the dominant backscatter mechanisms in nine inversions. 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. Section2 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.


## 2. 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.


### 2.1. 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.Figure 1
Definition of angular coordinates.The in-water wave vectors of incident and scattered directions are defined as follows:(1)ki=kwcos⁡θicos⁡∅iex+cos⁡θisin⁡∅iey-sin⁡θiez,(2)ks=kwcos⁡θscos⁡∅sex+cos⁡θssin⁡∅sey-sin⁡θies,where kw=2πf/cw is the acoustic wavenumber in water. cw and f are velocity of sound in water and acoustic frequency. Thus, the horizontal components and vertical component of the above wave vectors, respectively, are(3)Ki=kwcos⁡θicos⁡∅iex+cos⁡θisin⁡∅iey,(4)kiz=ki-Ki,(5)Ks=kwcos⁡θscos⁡∅sex+cos⁡θssin⁡∅sey,(6)ksz=ks-Ks,where the uppercase, boldface symbols represent the horizontal components. The following wave vector differences needed in the model calculation are defined:(7)ΔK=Ks-Ki,(8)Δkz=ksz-kiz.Based on Snell’s law, the wave vectors for the incident and scattered compressional waves in the sediment are(9)kpi=Ki-1ap2-cos2⁡θikwez,(10)kps=Ks-1ap2-cos2⁡θskwez,where ap 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(11)Δkp=Re⁡kps-kpi.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 ∅i=∅s=0, θs=π-θi for backscatter.


#### 2.1.1. 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 w2, ratio vp of sound velocity in sediment to cw, ratio aρ of sediment density to water density and, ratio δp of the imaginary and real parts of complex compressional speed. The final expression of scattering strength Sb and calculation process are as follows [17]:(12)Sb=10log10⁡kw4Aww2Ik2πΔK2Δkz2.Assuming that the roughness follows pure power law, the Kirchhoff integral Ik is(13)Ik=∫0∞J0ue-qu2αudu,where J0 is the zeroth-order cylindrical Bessel function of the first kind and α and q are as follows:(14)α=γ22-1,(15)q=12Ch2Δkz2ΔK-2α.Ch2 is the square of the “structure constant” which is related to the parameters of the power-law spectrum through(16)Ch2=2πw2Γ2-α2-2αα1-αΓ1+α,where Γ is the gamma function. The final expression of factor Aww is(17)Aww=121+Vwwθi1+VwwθsG,(18)G=1-1aρcos⁡θicos⁡θscos⁡∅s-1-ap2cos2⁡θi1-ap2cos2⁡θsap2aρ-1+1ap2aρ,where Vww is the flat-interface reflection coefficient,(19)Vwwθw=aρapsin⁡θw-1-ap2cos2⁡θwaρapsin⁡θw+1-ap2cos2⁡θw,where ap is the acoustic complex speed ratio which is defined as(20)ap=vp1+iδp.So far, the scattering strength Sb can be finally obtained.


#### 2.1.2. 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 withcw 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.vp, aρ, and δp have been defined before. The other three parameters are density fluctuations spectral exponent γ3, density fluctuations spectral strength w3, 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 Sv is [22, 23](21)Sv=10log10⁡1+Vwwθi21+Vwwθs2σv2kwap2Im1/ap2-cos2⁡θi+1/ap2-cos2⁡θs.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 σv through [24](22)σv=π2kw4μap2+cos⁡θicos⁡θscos⁡∅s-1ap2-cos2⁡θi1ap2-cos2⁡θs2WρρΔkp,where Wρρ is the spectrum for density fluctuations with the pure power law form:(23)Wρρk=w3kγ3.At last, the scattering strength of roughness and volume scattering (mixed scattering model) can be obtained through(24)Sall=10log10⁡1+Vwwθi21+Vwwθs2σv2kwap2Im1/ap2-cos2⁡θi+1/ap2-cos2⁡θs+kw4Aww2Ik2πΔK2Δkz2.


### 2.2. Bayesian Inference for Relative Entropy

A more detailed description of Bayesian inference applied to different inverse problems can be referred to [25–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) d and the scattering model are determined, the PPD of parameter vector m can be written as(25)Pm∣d=Pd∣mPmPd,where P(m) is the prior probability distribution of m and Pd∣m is conditional probability distribution of observations and is generally considered as likelihood function. P(d) is independent with m and is usually considered as constant. Then the PPD can be written as(26)Pm∣d∝Pd∣mPm=LmPm.Assuming zero-mean Gaussian-distributed errors, the likelihood can be expressed as(27)Lm=12πNd/2Cd1/2exp⁡-12d-dmTCd-1d-dm,where Nd is the dimension of d, d(m) denotes the vector predicted by the scattering model, and Cd is a diagonal data covariance matrix. One- and two-dimensional marginal PPD can be obtained as follows:(28)Pmi∣d=∫δmi′-miPm′∣ddm′,(29)Pmi,mj∣d=∫δmi′-miδmj′-mjPm′∣ddm′,where the non-boldface symbols mi, mj represent members of vector m. For the prior probability distribution and PPD, the expression of relative entropy is (30)Re=∫ln⁡Pm∣dPmPm∣ddm.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(31)I=∫fmPm∣ddm.Based on the theory of Monte Carlo, if the integral does not have an analytic solution, we can rewrite (31) as follows:(32)I=∫fmPm∣dQmQmdm=EQmfmPm∣dQm=EQmf′m.So the integral becomes the mean of f′(m) on the distribution function Qm. If we get a large number of samples m1,m2…mN from Qm, the integral can be approximately equal to(33)I≈1N∑i=1Nf′mi.If samples are generated from Pm∣d, the calculation process can be further simplified by(34)I≈1N∑i=1Nfmi.In some literatures, multiple maximum posterior (MAP) estimation methods [28, 29] are used to get samples with unknown Qm 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.


### 2.3. 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 statemc as the next state of mt-1 within Markov chain is(35)Prmt-1→mc=min⁡1,Pmcqmc→mt-1Pmt-1qmt-1→mc,where P(mc) and P(mt-1) are the posterior probability density of parameter vectors at different time and q(mc→mt-1) and q(mt-1→mc) are the conditional probabilities of chain moves from mc to mt-1 and mt-1 to mc, respectively. For DREAM, the conditional probability distribution is symmetric and the accepted probability of candidate state can be written as(36)Prmt-1→mc=min⁡1,PmcPmt-1.Then a random sample κ taken from a uniform distribution U(0,1) will be compared with Prmt-1→mc. The candidate state mc will be accepted as the next state of mt-1 if κ≤Pr. Otherwise, the next state mt will remain the same with mt-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.ii=1,2…,N chain is updated by(37)dmPi=εNp+1Np+λNpγδ,Np∑j=1δmPaj-mPbj,(38)dm≠Pi=0,where P is a subset (Np dimensional) of parameter space (Nm dimensional). Equation (37) in addition to (38) means there are only Np parameters changed by dmPi and keep the rest Nm-Np parameters unchanged within the No.i chain. The values of vector εNp and λNp are sampled independently from a multivariate normal distribution NNp(0,c) and a multivariate uniform distribution UNp(-c′,c′), respectively, with default values c′=0.1 and c=10-12. δ denotes the number of chains used to generate the perturbation. aj and bj denote the No.j element of vectors a and b consisting of δ integers drawn without replacement from 1,2,…,N. γδ,Np=2.38/2δNp 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. i chain at t is (39)mci=mt-1i+dmi.


## 2.1. 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.Figure 1
Definition of angular coordinates.The in-water wave vectors of incident and scattered directions are defined as follows:(1)ki=kwcos⁡θicos⁡∅iex+cos⁡θisin⁡∅iey-sin⁡θiez,(2)ks=kwcos⁡θscos⁡∅sex+cos⁡θssin⁡∅sey-sin⁡θies,where kw=2πf/cw is the acoustic wavenumber in water. cw and f are velocity of sound in water and acoustic frequency. Thus, the horizontal components and vertical component of the above wave vectors, respectively, are(3)Ki=kwcos⁡θicos⁡∅iex+cos⁡θisin⁡∅iey,(4)kiz=ki-Ki,(5)Ks=kwcos⁡θscos⁡∅sex+cos⁡θssin⁡∅sey,(6)ksz=ks-Ks,where the uppercase, boldface symbols represent the horizontal components. The following wave vector differences needed in the model calculation are defined:(7)ΔK=Ks-Ki,(8)Δkz=ksz-kiz.Based on Snell’s law, the wave vectors for the incident and scattered compressional waves in the sediment are(9)kpi=Ki-1ap2-cos2⁡θikwez,(10)kps=Ks-1ap2-cos2⁡θskwez,where ap 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(11)Δkp=Re⁡kps-kpi.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 ∅i=∅s=0, θs=π-θi for backscatter.


### 2.1.1. 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 w2, ratio vp of sound velocity in sediment to cw, ratio aρ of sediment density to water density and, ratio δp of the imaginary and real parts of complex compressional speed. The final expression of scattering strength Sb and calculation process are as follows [17]:(12)Sb=10log10⁡kw4Aww2Ik2πΔK2Δkz2.Assuming that the roughness follows pure power law, the Kirchhoff integral Ik is(13)Ik=∫0∞J0ue-qu2αudu,where J0 is the zeroth-order cylindrical Bessel function of the first kind and α and q are as follows:(14)α=γ22-1,(15)q=12Ch2Δkz2ΔK-2α.Ch2 is the square of the “structure constant” which is related to the parameters of the power-law spectrum through(16)Ch2=2πw2Γ2-α2-2αα1-αΓ1+α,where Γ is the gamma function. The final expression of factor Aww is(17)Aww=121+Vwwθi1+VwwθsG,(18)G=1-1aρcos⁡θicos⁡θscos⁡∅s-1-ap2cos2⁡θi1-ap2cos2⁡θsap2aρ-1+1ap2aρ,where Vww is the flat-interface reflection coefficient,(19)Vwwθw=aρapsin⁡θw-1-ap2cos2⁡θwaρapsin⁡θw+1-ap2cos2⁡θw,where ap is the acoustic complex speed ratio which is defined as(20)ap=vp1+iδp.So far, the scattering strength Sb can be finally obtained.


### 2.1.2. 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 withcw 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.vp, aρ, and δp have been defined before. The other three parameters are density fluctuations spectral exponent γ3, density fluctuations spectral strength w3, 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 Sv is [22, 23](21)Sv=10log10⁡1+Vwwθi21+Vwwθs2σv2kwap2Im1/ap2-cos2⁡θi+1/ap2-cos2⁡θs.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 σv through [24](22)σv=π2kw4μap2+cos⁡θicos⁡θscos⁡∅s-1ap2-cos2⁡θi1ap2-cos2⁡θs2WρρΔkp,where Wρρ is the spectrum for density fluctuations with the pure power law form:(23)Wρρk=w3kγ3.At last, the scattering strength of roughness and volume scattering (mixed scattering model) can be obtained through(24)Sall=10log10⁡1+Vwwθi21+Vwwθs2σv2kwap2Im1/ap2-cos2⁡θi+1/ap2-cos2⁡θs+kw4Aww2Ik2πΔK2Δkz2.


## 2.1.1. 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 w2, ratio vp of sound velocity in sediment to cw, ratio aρ of sediment density to water density and, ratio δp of the imaginary and real parts of complex compressional speed. The final expression of scattering strength Sb and calculation process are as follows [17]:(12)Sb=10log10⁡kw4Aww2Ik2πΔK2Δkz2.Assuming that the roughness follows pure power law, the Kirchhoff integral Ik is(13)Ik=∫0∞J0ue-qu2αudu,where J0 is the zeroth-order cylindrical Bessel function of the first kind and α and q are as follows:(14)α=γ22-1,(15)q=12Ch2Δkz2ΔK-2α.Ch2 is the square of the “structure constant” which is related to the parameters of the power-law spectrum through(16)Ch2=2πw2Γ2-α2-2αα1-αΓ1+α,where Γ is the gamma function. The final expression of factor Aww is(17)Aww=121+Vwwθi1+VwwθsG,(18)G=1-1aρcos⁡θicos⁡θscos⁡∅s-1-ap2cos2⁡θi1-ap2cos2⁡θsap2aρ-1+1ap2aρ,where Vww is the flat-interface reflection coefficient,(19)Vwwθw=aρapsin⁡θw-1-ap2cos2⁡θwaρapsin⁡θw+1-ap2cos2⁡θw,where ap is the acoustic complex speed ratio which is defined as(20)ap=vp1+iδp.So far, the scattering strength Sb can be finally obtained.


## 2.1.2. 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 withcw 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.vp, aρ, and δp have been defined before. The other three parameters are density fluctuations spectral exponent γ3, density fluctuations spectral strength w3, 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 Sv is [22, 23](21)Sv=10log10⁡1+Vwwθi21+Vwwθs2σv2kwap2Im1/ap2-cos2⁡θi+1/ap2-cos2⁡θs.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 σv through [24](22)σv=π2kw4μap2+cos⁡θicos⁡θscos⁡∅s-1ap2-cos2⁡θi1ap2-cos2⁡θs2WρρΔkp,where Wρρ is the spectrum for density fluctuations with the pure power law form:(23)Wρρk=w3kγ3.At last, the scattering strength of roughness and volume scattering (mixed scattering model) can be obtained through(24)Sall=10log10⁡1+Vwwθi21+Vwwθs2σv2kwap2Im1/ap2-cos2⁡θi+1/ap2-cos2⁡θs+kw4Aww2Ik2πΔK2Δkz2.


## 2.2. Bayesian Inference for Relative Entropy

A more detailed description of Bayesian inference applied to different inverse problems can be referred to [25–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) d and the scattering model are determined, the PPD of parameter vector m can be written as(25)Pm∣d=Pd∣mPmPd,where P(m) is the prior probability distribution of m and Pd∣m is conditional probability distribution of observations and is generally considered as likelihood function. P(d) is independent with m and is usually considered as constant. Then the PPD can be written as(26)Pm∣d∝Pd∣mPm=LmPm.Assuming zero-mean Gaussian-distributed errors, the likelihood can be expressed as(27)Lm=12πNd/2Cd1/2exp⁡-12d-dmTCd-1d-dm,where Nd is the dimension of d, d(m) denotes the vector predicted by the scattering model, and Cd is a diagonal data covariance matrix. One- and two-dimensional marginal PPD can be obtained as follows:(28)Pmi∣d=∫δmi′-miPm′∣ddm′,(29)Pmi,mj∣d=∫δmi′-miδmj′-mjPm′∣ddm′,where the non-boldface symbols mi, mj represent members of vector m. For the prior probability distribution and PPD, the expression of relative entropy is (30)Re=∫ln⁡Pm∣dPmPm∣ddm.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(31)I=∫fmPm∣ddm.Based on the theory of Monte Carlo, if the integral does not have an analytic solution, we can rewrite (31) as follows:(32)I=∫fmPm∣dQmQmdm=EQmfmPm∣dQm=EQmf′m.So the integral becomes the mean of f′(m) on the distribution function Qm. If we get a large number of samples m1,m2…mN from Qm, the integral can be approximately equal to(33)I≈1N∑i=1Nf′mi.If samples are generated from Pm∣d, the calculation process can be further simplified by(34)I≈1N∑i=1Nfmi.In some literatures, multiple maximum posterior (MAP) estimation methods [28, 29] are used to get samples with unknown Qm 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.


## 2.3. 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 statemc as the next state of mt-1 within Markov chain is(35)Prmt-1→mc=min⁡1,Pmcqmc→mt-1Pmt-1qmt-1→mc,where P(mc) and P(mt-1) are the posterior probability density of parameter vectors at different time and q(mc→mt-1) and q(mt-1→mc) are the conditional probabilities of chain moves from mc to mt-1 and mt-1 to mc, respectively. For DREAM, the conditional probability distribution is symmetric and the accepted probability of candidate state can be written as(36)Prmt-1→mc=min⁡1,PmcPmt-1.Then a random sample κ taken from a uniform distribution U(0,1) will be compared with Prmt-1→mc. The candidate state mc will be accepted as the next state of mt-1 if κ≤Pr. Otherwise, the next state mt will remain the same with mt-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.ii=1,2…,N chain is updated by(37)dmPi=εNp+1Np+λNpγδ,Np∑j=1δmPaj-mPbj,(38)dm≠Pi=0,where P is a subset (Np dimensional) of parameter space (Nm dimensional). Equation (37) in addition to (38) means there are only Np parameters changed by dmPi and keep the rest Nm-Np parameters unchanged within the No.i chain. The values of vector εNp and λNp are sampled independently from a multivariate normal distribution NNp(0,c) and a multivariate uniform distribution UNp(-c′,c′), respectively, with default values c′=0.1 and c=10-12. δ denotes the number of chains used to generate the perturbation. aj and bj denote the No.j element of vectors a and b consisting of δ integers drawn without replacement from 1,2,…,N. γδ,Np=2.38/2δNp 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. i chain at t is (39)mci=mt-1i+dmi.


## 3. 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 Section2.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 (40)d=dm+0.05×dm∘SD,where d(m) denotes the vector of true predicted scattering strength and d denotes the observation data vector, which is used for inversions. Each element of SD is drawn from the standard normal distribution. “∘” represents the operation: (A∘B)i,j=Ai,j(B)i,j. The angular and frequency ranges for d(m) are 10°~85° and 20kHz,40kHz,60kHz, respectively. Each of the three observation data sets is inverted once assuming each of the three types of scattering.Table 1
True parameter values for simulated inversions.
Parameter Parameter Parameter set_1 set_2 set_3 (roughness) (volume) (mixed) γ 2 3.25 3.25 3.25 w 2 0.000141 0.000141 0.000141 v p 1.1782 0.98 1.002 a ρ 2 1.39 1.66 δ p 0.01624 0.0011 0.0055 γ 3 3 3.63 3.35 w 3 0.0001 0.0137 0.00537 μ -1 1.667 1.389Figure2 shows the marginal PPD of each parameter from the inversion results for the true roughness scattering data. The red line, green line, and blue dotted line represent the inversion results based on roughness scattering model, volume scattering model, and mixed scattering model, respectively. It can be seen that inversion results of the mixed scattering model are better than that of the other two models on the whole. For the common parameters vp, aρ, and δp, inversion results of roughness scattering model are better than that of volume scattering model. The marginal PPDs of γ3, w3, and μ do not concentrate near the truth value; some even deviate from the truth value due to the dominant backscatter mechanism. Thus, we can probably judge that the dominant backscatter mechanism of this observation data is roughness scattering.Figure 2
Marginal PPD of each parameter assuming each of the three types of scattering (Parameter set_1).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 Figure3. It can be seen that the trend of uncertainty is consistent with the true scattering strength. Compared with large (θi>60°) and small (θi<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].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).The inversion results in Figure4 for true volume scattering are similar with the results in Figure 2. The 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, w2, which describe the interface roughness, are not properly estimated. The inversion result of aρ 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.Figure 4
Marginal PPD of each parameter assuming each of the three types of scattering (Parameter set_2).There is no drop for the predicted scattering strength of muddy seabed in Figure5 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 5
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_2).Figure6 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 w3 are relatively poor compared with the results in Figure 4.Figure 6
Marginal PPD of each parameter assuming each of the three types of scattering (Parameter set_3).Figure 7
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_3).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 Table2. 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 w3 in Figure 6 are relatively poor compared with the results in Figure 4.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.
roughness volume mixed (true) (true) (true) roughness -56.9702 -107.1008 -92.8925 volume -175.5403 -78.8034 -164.5447 mixed -53.9795 -78.8218 -76.5462


## 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 estimation. Three mechanisms for seabed scattering are considered: scattering from a rough water-seabed interface, scattering from volume heterogeneities, and mixed scattering from both interface roughness and volume heterogeneities. High-frequency seafloor roughness scattering and volume scattering are modelled based on Fluid Theories using Small-Slope Approximation and Small-Perturbation Fluid Approximation without considering the influence of the stratification 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.



---

*Source: 4272436-2018-12-02.xml*