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.
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. 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.
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.
The in-water wave vectors of incident and scattered directions are defined as follows:where 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, arewhere 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 arewhere is the acoustic complex speed ratio and its expression will be given in the next part. Then the in-sediment three-dimensional Bragg vector isIn 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 , 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 , roughness spectral strength , ratio 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 iswhere is the zeroth-order cylindrical Bessel function of the first kind and and are as follows: is the square of the “structure constant” which is related to the parameters of the power-law spectrum throughwhere is the gamma function. The final expression of factor iswhere is the flat-interface reflection coefficient,where is the acoustic complex speed ratio which is defined asSo far, the scattering strength 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 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. , , and have been defined before. The other three parameters are density fluctuations spectral exponent , density fluctuations spectral strength , 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 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 . Here the Small-Perturbation Fluid Approximation, which is also called Born approximation, is used to calculate through [24]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
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) and the scattering model are determined, the PPD of parameter vector can be written aswhere 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 asAssuming zero-mean Gaussian-distributed errors, the likelihood can be expressed aswhere 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 formBased 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 ,… from , the integral can be approximately equal toIf samples are generated from , the calculation process can be further simplified byIn 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.
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 state as the next state of within Markov chain iswhere and are the posterior probability density of parameter vectors at different time and and are the conditional probabilities of chain moves from to and to , respectively. For DREAM, the conditional probability distribution is symmetric and the accepted probability of candidate state can be written asThen a random sample taken from a uniform distribution will be compared with . The candidate state will be accepted as the next state of if . Otherwise, the next state will remain the same with . 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. chain is updated bywhere is a subset ( dimensional) of parameter space ( dimensional). Equation (37) in addition to (38) means there are only 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 and a multivariate uniform distribution , respectively, with default values and . 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 . is the jump rate which is periodically with a chance to enhance the probability of jumps between disconnected modes of the target distribution. Then the candidate state of No. chain at is
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 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 , 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 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 , , and , inversion results of roughness scattering model are better than that of volume scattering model. The marginal PPDs of , , 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.
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 () and small () 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 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. 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 , , 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 and 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 and in Figure 6 are relatively poor compared with the results in Figure 4.
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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was supported by the National Key Research and Development Program of China [No. 2016YFC1401203, No. 2018YFC1407400] and the National Natural Science Foundation of China [No. 41706106].