A Comparison of Three Sediment Acoustic Models Using Bayesian Inversion and Model Selection Techniques

Many geoacoustic models are used to establish the relationship between the physical and acoustic properties of sediments. In this work, Bayesian inversion and model selection techniques are applied to compare combinations of three geoacoustic models and corresponding scattering models—the fluid model with the effective density fluid model (EDFM), the grain-shearing elastic model with the viscosity grain-shearing (VGS(λ)) model, and the poroelastic model with the corrected and reparametrized extended Biot–Stoll (CREB) model. First, the resolution and correlation of parameters for the three models are compared based on estimates of the posterior probability distributions (PPDs), which are obtained by Bayesian inversion using the backscattering strength data. Then, model comparison and selection techniques are utilized to assess the matching degree of model predictions and measurements qualitatively and to ascertain the Bayes factors in favor of each quantitatively. These studies indicate that the fluid and poroelastic models outperform the grain-shearing elastic model, in terms of both parameter resolution and the ability to produce predictions in agreement with measurements for sandy sediments. The poroelastic model is considered to be the best, as the inversion based on it can provide more highly resolved information of sandy sediments. Finally, the attempt to implement geoacoustic inversion with different models provides a relatively feasible remote sensing scheme for various types of sediments under unknown conditions, which needs further validation.


Introduction
As the main means of underwater investigation, acoustic remote sensing techniques have been a research hotspot over the last decades.Underwater topographic acoustic telemetry has been a great success, leading to the emergence of many highly efficient commercial products, ranging from high-frequency single-and multibeam echosounders and sidescan sonars to low-frequency sub-bottom profiling tools and seismic systems.In the course of investigation using such equipment, sediment sampling has often been done to form a dataset corresponding to the acoustic data, and the sediment samples are used for the verification of seabed classification [1].Restricted by the time and cost required for collection and analysis of samples, the verification behavior cannot be expanded on a large scale.Indeed, seabed classification just tells us that the sediments are different at different locations.
In order to obtain geoacoustic properties of sediments, model-based inversion of acoustic data is a more effective and feasible means compared with grab samples [2].
According to the sediment geoacoustic model categories proposed by Jackson and Richardson [3], there are three categories, which are fluid, elastic, and poroelastic, based on the number and types of acoustic waves that can travel in sediments.We can see from the literature [4][5][6][7][8][9] that these geoacoustic models are mainly used to establish the relationship between physical and acoustic properties of sediments, such as sound speed and attenuation, and are regarded as "wave theory" or "propagation theory" to distinguish them from subsequent scattering models by Jackson and Richardson.In acoustic remote sensing, scattering models of water-seabed interface are necessary to calculate observations such as scattering strength in model-based inversion, since the speed and attenuation of sound in sediments cannot be directly measured from a distance.Similarly, scattering models are also classified into three categories.
Due to the prevalence of sandy sediment in coastal environments, there have been a relatively large number of experimental data/model comparisons corresponding to sandy sediment [10][11][12][13].The development of models is driven by the mismatch of model prediction and acoustic observation, including sound speed and attenuation, within sediments measured in situ and scattering strength measured remotely.The poroelastic model may be a reasonable choice for a porous medium such as sandy sediment.Biot first proposed a theory of sound propagation in a poroelastic medium, taking into account the interactions between the pore fluid and consolidated elastic frame [14,15].The theory was first used to model sandy sediments by Stoll [16].The Biot-Stoll model is a fundamental improvement over the fluid and viscoelastic models and provides a satisfactory fitting performance of model prediction with the measured data [10,13].Based on the mismatch between model prediction and measured data, Chotiros further extended the Biot-Stoll model to include the physics of grain-grain contact, multiple scattering losses, a more efficient set of input parameters, and a heuristic correction that allows the model to span a wide range of sediment types [17], which was redefined as the corrected and reparametrized extended Biot-Stoll (CREB) model by Bonomo and Isakson [18].
In addition to direct modification of the Biot-Stoll model, Williams presented an acoustic propagation model that approximates a porous medium as a fluid with a bulk modulus and effective density derived from Biot-Stoll theory [19].It was shown that the dispersion, transmission, reflection, and scattering predicted with this approximate model, which is named the effective density fluid model (EDFM), are very close to the predictions of the full Biot-Stoll theory for sandy sediment [10,12,20].The approximate effectiveness and limits of EDFM were analyzed by Bonomo and co-workers using the finite element method (FEM) [21].In order to further improve the matching degree between the measurements and the model predictions [20], two additional physical mechanisms, transfer of heat between the liquid and solid and the effect of granularity, were added to the EDFM [7].
Although the poroelastic model and its approximation are successful to some extent, another kind of geoacoustic model has received considerable attention.Based on the assumption that the shear rigidity modulus of the medium is zero and implying the absence of a skeletal elastic frame, Buckingham proposed the grain-shearing (GS) elastic model [22].The model was extended to include the effects of viscosity of the molecularly thin layer of pore fluid separating contiguous grains [6].The extended version is named VGS, with the first initial serving as a reminder that the model includes viscosity of the pore fluid.The latest version of Buckingham's model is designated the VGS(λ), with a modification to ensure that the viscosity of the pore fluid has no effect on the shear wave.The shear dispersion relations in the VGS and GS theories are identical.At the same time, Buckingham showed the adequate matching performance of the model prediction with the measured data [8,23].
The main content of the current paper is to compare the performance of inversions based on the three models, EDFM, VGS(λ), and CREB, and their corresponding acoustic scattering models.Bayesian model selection techniques are used to evaluate the degree of matching between model predictions and measurements, which are backscattering strengths in this paper.In fact, this is not the first time that Bayesian techniques have been applied to underwater acoustics.Dosso and co-workers successfully applied Bayesian techniques to estimate the number of seabed layers and the uncertainty of seabed parameters [24][25][26].There have also been attempts to use Bayesian techniques to solve the problem of seabed sediment classification [27,28].Bayesian inversion and model selection based on the geoacoustic models using sound speed and attenuation measurements were carried out by Bonomo and Isakson to compare the competing geoacoustic models [18].In this work, we introduce a similar comparison for underwater remote sensing by combining the geoacoustic models with corresponding scattering models, using scattering strength measurements to carry out inversion and model selection.The authors believe this is the first time that model selection techniques have been applied to underwater remote sensing using scattering strength to compare different types of models.Before parameter inversion, a parameter sensitivity analysis is carried out to forecast the relative resolution of parameters of each model and provide an interpretation of the subsequent inversion results.Then, the resolution of parameters that are mostly difficult to measure in situ can be obtained from inversion results and compared with the results of the sensitivity analysis.The model selection techniques can provide a quantitative measure of the most suitable model.When faced with a complex and unknown sediment environment, the comparison can also provide a relatively feasible scheme for assessing different kinds of models for underwater acoustic inversion and remote sensing, and the best model and inversion results for specific sediment can be obtained.
The observation used in the inversion is backscattering strength varying with angle and frequency.Since this paper is mainly concerned with the uppermost layer of sediment and the scattering model used here is only suitable for high frequency, the backscattering strength measurements are carried out at high frequencies (≥100 kHz).The remainder of this paper has the following organization.The calculation processes of scattering strength for the three models and corresponding model parameter sensitivity analysis are given in Section 2. Section 3 briefly reviews the Bayesian inference method used in this work.Section 4 presents the experimental measurements, parameter inversion, and model selection results.Finally, this work is concluded in Section 5.

Geoacoustic Models and Acoustic Scattering Models
Geoacoustic models are usually established to describe the relationship between physical parameters and the attenuation and speed of sound waves that propagate within sediment, as mentioned above.The main content of this section is a combination of the latest geoacoustic and scattering models used to calculate scattering strength.According to the model categories proposed by Jackson and Richardson [3], a fluid medium only supports the propagation of compressional waves, an elastic medium supports the propagation of compressional waves and shear waves, and a poroelastic medium supports additional slow compressional waves compared to the elastic medium.It is important to note that slow compressional waves are given theoretically and have never been observed in natural marine sediment.The combination means that the attenuation and speed of sound waves supported by the corresponding geoacoustic model categories are needed as input for different kinds of scattering models.Indeed, different kinds of geoacoustic and scattering models can be applied in cross-combinations to calculate scattering strength.The EFDM is an example in which effective density derived from the Biot-Stoll theory is used to calculate scattering strength based on fluid scattering models.Then the fluid scattering model is extended to the poroelastic case by the simple replacement of density by effective density.
Although the separation of scattering into separate, independent roughness and volume components is somewhat artificial, it is helpful to understand the scattering mechanisms by discussing them separately.The total scattering strength can be expressed as where σ bI and σ bV denote scattering cross-sections of rough water-seabed interface and volume heterogeneities, respectively.As mentioned above, due to the prevalence of sandy sediment, the study in this paper is also based on sandy sediment.At high frequencies, the scattering of sandy sediment is mainly due to interface roughness scattering [10], which is modeled by three types of models, and the contribution from a fluid volume scattering model is included for completeness in this paper.Through the cross-combination mentioned above, the fluid volume scattering model can also be extended to the elastic and poroelastic case.The following are the calculation processes of scattering strength of the fluid model, grain-shearing elastic model, and poroelastic model, and the corresponding model parameter sensitivity analysis.

Fluid Model
Although seabed sediments are composed of discrete particles, it is a reasonable approximation that they are regarded as a continuous medium [29] and this assumption applies to all models used in this paper.As a relatively simple approach, fluid theories are commonly used to describe sediment acoustic properties and have been combined with scattering models.The EFDM, which is proposed as a fluid approximation of the Biot-Stoll theory, extends the fluid model input parameters to the poroelastic case.The calculation process of this approximation only considers the propagation of compressional waves in sediments, leading to its category as a fluid type.The fluid model used here to calculate the scattering strength requires 15 input parameters-roughness spectral exponent γ 2 , roughness spectral strength ω 2 , density fluctuation spectral exponent γ 3 , density fluctuation spectral strength ω 3 , ratio of compressibility to density fluctuation in sediment µ, mean grain diameter d, tortuosity α, porosity β, dynamic viscosity of water η, permeability κ, mass density of water ρ w , mass density of grains ρ g , bulk modulus of water K w , bulk modulus of grains K g , and compressional wave speed in water c w .

EDFM
The fluid model used in this paper has been extended to the poroelastic case by using the effective density.Therefore, the calculation of scattering strength needs to begin from the relevant parameters of the poroelastic theory.Based on the Biot constitutive relations, four moduli are given successively by Stoll, and the relationships between the moduli and the parameters of constituent media are as follows [30]: (2) where K f is the complex frame bulk modulus and U is the complex frame shear moduli.The starting point of the EDFM regarded as an approximation for full Biot-Stoll is the frame modulus set to zero.The reason for this is that the frame moduli are small and can be negligible compared with the grain and water moduli.Then an effective fluid modulus can be obtained: The effective fluid modulus is substituted into the wave equation of the Biot theory, and the detailed derivation process can be found in [19].The effective density, complex acoustic compressional speed, and corresponding phase speed in sediment can be obtained by: where ρ = αρ w /β + iFη/(2π f κ), ρ = βρ w + (1 − β)ρ g , f is the acoustic frequency, and F is the complex correction factor for dynamic viscosity: where J 0 and J 1 are Bessel functions of the first kind, and a = 8ακ/β is called the pore size.In order to further improve the matching degree between model predictions and measurements of sound speed and attenuation, Williams added two additional nongrain contact mechanisms to the EDFM: heat transfer between the liquid and solid at low frequencies and granularity of the medium at high frequencies [7].Since the application of the EDFM in this paper is at high frequencies, only the model modification using doublet mechanics (DM) to model the granularity for high frequency is given here.This is done by adding a correction factor ξ DM to K e f f as follows: where λ = c w / f is the wavelength of sound wave in water and ξ DM can be obtained by: where T i1 are direction cosines and functions of θ, and the expressions are: The additional multiple scattering loss α ms = 24(2π f d/c w ) 4 is added into the attenuation (dB/m) by [31]: The final modified complex compressional wave speed can be expressed as: The complex speed ratio a 1 and density ratio a ρ needed to calculate the scattering strength in the scattering model are:

Fluid Interface Roughness Scattering Model
For roughness scattering at the sea bottom interface, two widely used approximations are used to calculate the scattering strength: the small-roughness perturbation method and the Kirchhoff approximation.The perturbation theory tends to be most accurate for scattering at wide angles relative to the specular (flat-interface reflection) direction [32], and the Kirchhoff approximation is better for scattering near the specular direction [33].The approximation method used in this paper is called small-slope approximation; it was first proposed by Voronovich to study scattering by the sea surface [34] and was subsequently used to study seabed scattering.In the interest of simplicity, the lowest-order small-slope approximation is used, 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 approximation.According to the summary by Jackson, the uniform expression of roughness scattering cross-section for different models using lowest-order small-slope approximation is [3]: where k w = 2π f /c w is the acoustic wavenumber in water and ∆K, ∆k z denote the magnitude of the horizontal and vertical vector difference between the scattered wave and the incident wave, respectively, and are given in Appendix A.
Assuming that the roughness is isotropic and follows pure power law, I k can be obtained by: where C 2 h 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 above calculation process of roughness scattering cross-section is also applicable to the grain-shearing elastic and poroelastic models.The difference is the calculation of A ww .The expression of factor A ww for the fluid model is: where θ i , θ s , and ∅ s are defined in Figure A1.V ww is the flat-interface reflection coefficient and can be obtained by: Z(θ) = a ρ a 1 sin θ/ 1 − a 2 1 cos 2 θ (25)

Fluid Volume Scattering Model
Although the acoustic scattering caused by the inhomogeneity of sediments is more significant only in soft seabeds, such as muddy seabed, where sound waves can easily travel in and out, the volume scattering is taken as a supplement and modification for total scattering strength to further enhance the matching performance of the models.The small-perturbation fluid approximation for volume scattering is used.Assuming that the volume scattering strength is independent of depth, the expression for σ bV is [35,36]: where ∆k p is the real part of the vector difference of the scattered and incident waves that propagate in sediment, given in Appendix A. So far, the scattering strength for the fluid model can be obtained through Equations ( 1), (18), and (26).

Grain-Shearing Elastic Model
The categories of models discussed by Jackson and Richardson are meant to be broad and can apply to consolidated sediments such as rock where model proposed by Buckingham is designed for unconsolidated sediments such as sand.In order to distinguish the elastic model used in this paper from the empirical visco-elastic model developed by Hamilton [37], the term 'grain-shearing elastic model' is used to specify the combination of VGS(λ) and elastic scattering model.Compared with the fluid model, the grain-shearing elastic model incorporates shear forces.The sediment supports the propagation of shear waves, which is modeled to provide additional corrections.The grain-shearing elastic model used here to calculate the scattering strength requires 16 input parameters: roughness spectral exponent γ 2 , roughness spectral strength ω 2 , density fluctuation spectral exponent γ 3 , density fluctuation spectral strength ω 3 , ratio of compressibility to density fluctuation in sediment µ, porosity β, mass density of water ρ w , mass density of grains ρ g , bulk modulus of water K w , bulk modulus of grains K g , material exponent n, compressional rigidity coefficient γ 1 , shear rigidity coefficient γ t , compressional viscoelastic relaxation time τ 1 , shear viscoelastic relaxation time τ t , and compressional wave speed in water c w .

VGS(λ)
As the latest version of the model proposed by Buckingham, VGS(λ) is an extension of GS that includes the effects of viscosity of the molecularly thin layer of pore fluid separating contiguous grains and a slight modification of the shear viscoelastic relaxation time.The VGS(λ) dispersion curves have also been verified by the measurements in SAX99 [8].
First, Wood's equation, which usually gives a rather poor prediction of sediment sound speed [38], is given, as it ignores much of the dynamics of the grain-pore water system: From the view of expression, the VGS(λ) expressions can be considered as modifications of c 0 and the compressional wave speed c 1phase , attenuation α 1 , shear wave speed c tphase , and attenuation α t , respectively, are expressed as [6]: where τ t is the shear viscoelastic relaxation time and τ t → ∞ is taken because of the experimental evidence [8] that τ t τ 1 .Then the complex compressional wave speed c 1 and shear wave speed c t can be similarly obtained by Equation (15).The complex speed ratio a 1 and a t can be obtained by Equation (16).The density ratio a ρ for the grain-shearing elastic model is as follows: Remote Sens. 2019, 11, 562 8 of 30

Elastic Scattering Model
The expression of scattering cross-section for the elastic interface roughness scattering model is similar to the one used with the fluid model.The difference is the factor A ww used in Equation (18), which can be obtained by [3]: where D i , i = 1, 2, 3, 4, given in Appendix B; V ww (θ) is the flat-interface reflection coefficient, which can be obtained by Equation (24); and the corresponding Z(θ) for the elastic scattering model is: As mentioned above, the fluid volume model can be extended for use in an elastic medium by substituting the appropriate parameters into Equations ( 26) and (27).Then the total scattering strength for the grain-shearing elastic model can be obtained by Equation (1).

Poroelastic Model
Both elasticity and porosity of sediments are considered in the poroelastic model, where the sediment grains constitute an elastic "frame" coupled to the pore fluid.This is the reason that the poroelastic model supports additional "slow" compressional waves compared with the elastic model.This model is essentially the most appropriate for sandy sediment, which is composed of grains and water as a two-phase system [16].The total poroelastic model requires 19 input parameters: roughness spectral exponent γ 2 , roughness spectral strength ω 2 , density fluctuation spectral exponent γ 3 , density fluctuation spectral strength ω 3 , ratio of compressibility to density fluctuation in sediment µ, mean grain diameter d, porosity β, dynamic viscosity of water η, mass density of water ρ w , mass density of grains ρ g , bulk modulus of water K w , bulk modulus of grains K g , cementation exponent m, pore shape factor a B , Poisson's ratio of grains ν, low-frequency asymptotic frame bulk modulus K bo , high-frequency asymptotic increase K y , bulk relaxation frequency f k , and compressional wave speed in water c w .

CREB
In order to reproduce the frequency dependence of the measured sound speed and attenuation and achieve a better matching degree between model predictions and measurements, grain contact physics models have been developed and tested for the poroelastic model.Chotiros and Isakson used grain contact squirt flow and viscous drag mechanisms to better model the observed behavior of attenuation [4].Then a missing relaxation mechanism was added at higher frequencies related to the grain contact physics, and the extended Biot model was proposed to achieve further improvement [9].The model used here is the latest version in which the extended Biot model is reparametrized and corrected to be able to model sediments of multiple types [17].
Permeability κ and tortuosity α can be related to the CREB parameters through: Based on the extended Biot model, the complex effective frame bulk modulus K f and frame shear modulus U are replaced by frequency-dependent expressions: where F g can be set equal to 1 for sandy sediment or to other values to model different types of sediments.The moduli H, C, M, complex correction factor F can be obtained by Equations ( 2)-( 4) and (10) and the "fast" and "slow" compressional wavenumbers can be expressed as the roots of the following expression: where root k 1 with a negative imaginary component and a smaller real component is the "fast" compressional wavenumber, and root k 2 is the "slow" compressional wavenumber.The shear wavenumber k t is the root with a negative imaginary component of the following expression: The density ratio a ρ for the poroelastic model is the same as that of the grain-shearing elastic model, and the speed ratio can be obtained by:

.2. Poroelastic Scattering Model
Similar to the elastic scattering model, the scattering cross-section for the poroelastic interface roughness scattering model can be obtained by Equation (18), and the corresponding factor A ww is the first element of the following column matrix [39]: which can be computed by: where the subscripts of Y are used to distinguish the two matrices without specific meaning and the superscripts of matrices E, P, and Q correspond to the three spatial components.The matrices Y, E, P, and Q are given in Appendix B. V 0 (K i ) is a six-row column vector comprising five reflection and transmission coefficients and supplemented with unity in the last row.The flat-interface reflection coefficient V ww (θ) for the poroelastic model is the first element of the column vector V 0 (K i ).Then the volume scattering can also be computed by Equations ( 26) and ( 27) and the total scattering strength for the poroelastic model can be obtained by Equation (1).

Sensitivity Analysis
Parameter sensitivity analysis is from the view of mathematics and is a global analysis method.The potential correlation between some parameters is not specified in advance and is not considered in the analysis process.Since model predictions of some parameter combinations may not be consistent with local measurements, the parameter resolution given by the results of the sensitivity analysis may not be exactly consistent with that given by inversion results.However, we can still obtain a preliminary understanding of the resolutions of the model parameters through sensitivity analysis.The sensitivity analysis can forecast the relative resolution of parameters for each model and provides an interpretation of the subsequent inversion results.Similar to many global sensitivity analysis methods, the uncertainty importance measure proposed by Borgonovo is used here [40].Borgonovo verified that the parameter that influences variance the most is not necessarily the parameter that influences the output distribution the most, and proved the effectiveness of his method over other variance-based methods [41][42][43][44].The definition of the first-order Borgonovo index is: where f Y is the probability distribution of the model prediction and f Y|X i is the conditional distribution on X i .δ i indeed is a measure of the expected shift in the probability distribution of the model prediction when X i is set to a fixed value.The parameter space is determined for sandy sediments based on the results from previous research in the literature.Bounded, uniform distributions are taken as the parameter priors, and the distribution parameters are given in Table 1, where they are divided into the common parameters of the three models and the unique parameters of each model.The priors are chosen to be minimally informative and represent attempts at capturing the full range of possible values the model input parameters can take for sandy sediments.As shown in Table 2, the parameters that describe the characteristics of water and the sound speed in water are supposed to be known parameters, which refer to water parameters listed in [45].The reason is that the subsequent experiments are carried out in a rectangular tank filled with fresh water.It should be noted that the mass density and bulk modulus of grains are replaced by ratios to water for convenience, as shown in Table 1.The sensitivity analysis here and the subsequent parameter inversions are both based on these priors.Using Latin hypercube sampling [46], 60,000 samples are generated to calculate the Borgonovo index for each model.The grazing angle range and frequency, respectively, are 20 • -70 • and 100 kHz.
Figure 1a displays the Borgonovo indices for the fluid model.Inspection of the figure reveals that the roughness spectral exponent γ 2 , density fluctuation spectral exponent γ 3 , and porosity β have a great influence on model prediction.The effect on the model prediction of roughness spectral strength ω 2 and ratio of compressibility to density fluctuations µ is smaller than that of the first three parameters but larger than that of all the remaining parameters.This means that there could be five parameters that yield better resolution in the fluid model-based inversion.Figure 1b shows that the roughness spectral exponent γ 2 and density fluctuation spectral exponent γ 3 still have a significant impact on model prediction for the grain-shearing elastic model, while roughness spectral strength ω 2 , ratio of compressibility to density fluctuation µ, porosity β, and material exponent n only have a moderate effect on model prediction.For the poroelastic model in Figure 1c, the sensitive parameters are roughly the same as the fluid model, and the others remain insensitive to model prediction.From Figure 1, we can also see that the indices of some parameters fluctuate under different grazing angles.Indeed, sensitivity analysis was also performed at different frequencies, and the results show that the frequency change had only a slight effect on the indices.Thus, we only give results of sensitivity analysis at 100 kHz here.parameters are roughly the same as the fluid model, and the others remain insensitive to model prediction.From Figure 1, we can also see that the indices of some parameters fluctuate under different grazing angles.Indeed, sensitivity analysis was also performed at different frequencies, and the results show that the frequency change had only a slight effect on the indices.Thus, we only give results of sensitivity analysis at 100 kHz here.

Bayesian Inference
Bayesian inference is by far the most suitable solution to the inversion problem, as it treats all the parameters as random variables.The inversion results are given in the form of posterior probability distribution (PPD) with quantized uncertainty instead of point estimation based on the global optimization algorithm.This section only makes a brief review; more complete treatments of Bayesian inference applied to geoacoustic inversion can be found elsewhere [24,25,[47][48][49].When a set of observations (scattering strength)  and the adopted model  for inversion are determined, the PPD of parameter vector  can be expressed as where () is the prior probability distribution of .(|, ) is the conditional probability distribution of observations  given model parameters  and model .() is independent of  and is usually considered as a constant factor.Then the PPD can be written as

Bayesian Inference
Bayesian inference is by far the most suitable solution to the inversion problem, as it treats all the parameters as random variables.The inversion results are given in the form of posterior probability distribution (PPD) with quantized uncertainty instead of point estimation based on the global optimization algorithm.This section only makes a brief review; more complete treatments of Bayesian inference applied to geoacoustic inversion can be found elsewhere [24,25,[47][48][49].When a set of observations (scattering strength) d and the adopted model Mod for inversion are determined, the PPD of parameter vector m can be expressed as P(m|d, Mod ) = P(d|m, Mod )P(m)/P(d) (48) where P(m) is the prior probability distribution of m.P(d|m, Mod ) is the conditional probability distribution of observations d given model parameters m and model Mod.P(d) is independent of m and is usually considered as a constant factor.Then the PPD can be written as P(m|d, Mod ) ∝ P(d|m, Mod )P(m) = L(m|Mod )P(m) (49) where P(d|m, Mod ) is also known as the likelihood function L(m|Mod ).In order to define the likelihood function, the error of observations, which consists of measurement error and model error, needs to be specified.An effective approach is used here in which the errors of observations are taken as unknown random parameters, like input parameters for model prediction, and estimated simultaneously within the inversion process.Assuming common independent Gaussian-distributed errors, the likelihood function can be expressed as where N d is the dimension of d; σ i is the standard deviation for d i , which is the ith element of d; and d i (m|Mod ) is the ith element of d(m|Mod ), which denotes the observations predicted by the parameter vector m and corresponding model Mod.The solution of the inversion problem is given in the form of a marginal PPD, defined as: where m i and m j represent the elements of parameter vector m and P(m i |d, Mod) and P m i , m j d, Mod denote the one-and two-dimensional marginal PPD, respectively.The integral that needs to be solved can be rewritten as: For nonlinear problems, such as geoacoustic inversion, the integral does not have an analytic solution.Based on Monte Carlo theory, we can rewrite the integral as: If a large number of samples m 1 , m 2 . . .m N s is obtained from Q(m), the integral can be approximately equal to: Let Q(m) = P(m|d, Mod ), then the integral can be further simplified to: which means that the main objective of inversion is to obtain a set of samples from the PPD to approximate the target integral.It is not difficult to figure out that a larger number of samples that follow the target probability distribution can better approximate the true value of the integral.The rigorous sampling methods will be given in the next section.

Parameter Inversion
Monte Carlo theory provides the idea that samples can be used to approximate an integral without an analytic expression.Then the key task of parameter inversion is to generate a set of samples that can fully represent the target probability distribution.Compared with the static simulation of traditional Monte Carlo integral, the Markov chain Monte Carlo (MCMC) method is a dynamic simulation as the process goes on.MCMC exploits a Markov chain that generates a random walk through the parameter search space to obtain solutions with stable frequencies stemming from a stationary distribution.Different MCMC methods have different constructions of transfer probability, which has a crucial influence on the convergence efficiency of the Markov chain.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 known as the Metropolis-Hastings (M-H) algorithm [50].In order to increase the sampling efficiency of complex posterior distributions involving long tails, correlated parameters, multimodality, and numerous local optima, more improved algorithms have been proposed, such as adaptive proposal (AP) [51], adaptive Metropolis (AM) [52], delayed rejection adaptive Metropolis (DRAM) [53], differential evolution Markov chain (DE-MC) [54], fast Gibbs sampling (FGS) [47], and so on.The sampling method used in this paper is a relatively efficient multichain method called differential evolution adaptive Metropolis (DREAM), which can maintain detailed balance and ergodicity with excellent performance on a wide range of problems involving nonlinearity, high dimensionality, and multimodality [55].In general, the probability of accepting candidate state m c as the next state of m t−1 within the Markov chain is: where P(m c ) and P(m t−1 ) are the posterior probability densities of parameter vectors at different times, and q(m c → m t−1 ) and q(m t−1 → m c ) are the conditional probabilities of chain moves from m c to m t−1 and m t−1 to m c , respectively.For DREAM, the conditional probability distribution is symmetric and the probability of accepting the candidate state can be written as: Then, a random sample S r taken from a uniform distribution U(0, 1) will be compared with P r (m t−1 → m c ).The candidate state m c will be accepted as the next state of m t−1 if S r ≤ P r .Otherwise, the next state m t will remain the same as m t−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 ith (i = 1, 2 . . ., N) chain is updated by: where P is a subset (N p dimensional) of parameter space (N m dimensional).Equations ( 59) and (60) mean there are only N p parameters changed by dm i P and the other N m − N p parameters remain unchanged within the ith chain.The values of vectors ε N p and λ N p are sampled independently from a multivariate normal distribution N N p (0, c) and a multivariate uniform distribution U N p (−c , c ), respectively, with default values c = 0.1 and c = 10 −12 .δ denotes the number of chains used to generate the perturbation.a j and b j denote the jth element of vector a and b consisting of δ integers drawn without replacement from {1, 2, . . . ,N}. γ (δ,N p ) = 2.38/ 2 δN p is the jump rate, which is periodical with a 20% chance to enhance the probability of jumps between disconnected modes of the target distribution.

Convergence Criterion
An important aspect when applying the MCMC method is determining when to safely stop.The solution proposed by Dosso [47,56] is that samples of two independent chains are periodically compared until they are deemed to have satisfied an objective convergence criterion.In this paper, the DREAM algorithm is adopted, which uses multiple chains to achieve different trajectories running in parallel and explore the posterior target distribution.Compared with the dual chain scheme used by Dosso, multiple chains are more applicable for complex posterior distributions involving long tails, correlated parameters, multimodality, and numerous local optima, offering robust protection against premature convergence and opening up the use of a wide arsenal of statistical measures to test whether convergence to a limiting distribution has been achieved [57].Here the convergence of the inversion process is monitored with the R statistic proposed by Gelman and Rubin [58,59].If R is large, this means the parameter space is not fully explored by the chains, and more iterative simulations are needed.Alternatively, if R is close to 1, it can be concluded that the samples within the chains are close to the target probability distribution.The threshold of R used in this paper is 1.2, and the actual convergence criterion is that R remains below the threshold for 3000 simulations.

Model Selection
One of the main purposes of this paper is to determine which model best matches the measured data.For model-based Bayesian inversion, the Bayes factor [60] is the most common evaluation criterion, which provides a quantitative index to judge competing models.The Bayes factor of Mod 1 with respect to the alternative Mod 2 can be calculated by: (61) which is equivalent to the ratio of the evidence, EV 1 and EV 2 , for two competing models.The evidence for a selected model Mod i is defined as the integral of the likelihood function over the prior distribution: In practice, most regions of the parameter space in the prior distribution correspond to small likelihood function values and make negligible contributions to the integral.Similar to importance sampling, the prior distribution can be replaced by the PPD to restrict the integral space to those areas that make significant contributions to the integral, reducing unnecessary parameter space sampling in abundance and computation to approximate the evidence to the greatest extent [61].Similar to Equations ( 53)-( 56), the samples within the Markov chains are used to approximatively calculate the evidence in this paper by: The value of evidence can be interpreted by the scale proposed by Robert and Adrian [60], as shown in Table 3, where the scale is divided into four (increasing) levels of support for proposition Mod 1 relative to Mod 2 .

Experimental Measurements, Results, and Discussion
This section presents and analyzes the results of the Bayesian inference performed in this work.First, a description of the laboratory measurements is given in Section 4.1.Then, the parameter inversion results based on the laboratory measurements are analyzed.Finally, the model selection and comparison results are discussed in Section 4.3.

Experiment Description
Specialized to the monostatic case, backscattering strength was measured in a rectangular tank using a directional transducer.As shown in Figure 2, a circular arc was used to mount the transducer, and the bottom was filled with river sand.The thickness of the sand layer in the tank was 0.5 m, and the diameter of the circular arc was 3 m.The associated data-processing methods have been described in detail by Jackson [3,62], and only a brief summary of the experimental configuration will be given here.The mean-square voltage at the terminals of the transducer was averaged over groups of 50 successive pings consisting of 10 measurements at five evenly spaced positions in the y direction as the arc moved along the y-axis.This double-averaging was done to mitigate the fluctuation of a single measurement and the potential weak heterogeneity differences of the sand as well as the unevenness of the bottom surface.The transducer was high-frequency broadband, which could move along the arc, and separate passes were made with different grazing angles and frequencies.The grazing angle range and frequencies used in the experiment were 20

Experiment Description
Specialized to the monostatic case, backscattering strength was measured in a rectangular tank using a directional transducer.As shown in Figure 2, a circular arc was used to mount the transducer, and the bottom was filled with river sand.The thickness of the sand layer in the tank was 0.5 m, and the diameter of the circular arc was 3 m.The associated data-processing methods have been described in detail by Jackson [3,62], and only a brief summary of the experimental configuration will be given here.The mean-square voltage at the terminals of the transducer was averaged over groups of 50 successive pings consisting of 10 measurements at five evenly spaced positions in the y direction as the arc moved along the y-axis.This double-averaging was done to mitigate the fluctuation of a single measurement and the potential weak heterogeneity differences of the sand as well as the unevenness of the bottom surface.The transducer was high-frequency broadband, which could move along the arc, and separate passes were made with different grazing angles and frequencies.The grazing angle range and frequencies used in the experiment were 20°-70° and 100 kHz, 120 kHz, and 140 kHz, respectively.There is a slight drop of backscattering strength above the critical angle of 30° which is consistent with the measurements given by Williams [10].Although no field experiments were carried out, the authors believe that laboratory measurements are more manageable and provide adequate fitness between model and data.Of course, more subsequent field data for model selection analysis are required due to the limitations of laboratory measurements.This work features the use of a remote sensing scheme with multiple models and lays a foundation for future practical field applications.There is a slight drop of backscattering strength above the critical angle of 30 • which is consistent with the measurements given by Williams [10].Although no field experiments were carried out, the authors believe that laboratory measurements are more manageable and provide adequate fitness between model and data.Of course, more subsequent field data for model selection analysis are required due to the limitations of laboratory measurements.This work features the use of a remote sensing scheme with multiple models and lays a foundation for future practical field applications.
frequencies.It can be seen that the backscattering strength has only a weak frequency dependence.There is a slight drop of backscattering strength above the critical angle of 30° which is consistent with the measurements given by Williams [10].Although no field experiments were carried out, the authors believe that laboratory measurements are more manageable and provide adequate fitness between model and data.Of course, more subsequent field data for model selection analysis are required due to the limitations of laboratory measurements.This work features the use of a remote sensing scheme with multiple models and lays a foundation for future practical field applications.

Parameter Inversion
The parameter priors used for the parameter inversion were the same as those for the sensitivity analysis as shown in Tables 1 and 2. A total of five chains were used during the sampling process and deemed sufficient for parameter dimensions of all models.Under the premise of final convergence, the number of total generations varied between 120,000 and 400,000 depending on the model.The PPDs and Bayes factors were estimated by Equations ( 51), (52), and (61) using the last 60,000 DREAM samples.Estimates of the one-and two-dimensional marginal PPDs for the three models obtained via Bayesian parameter inversion applied to the backscattering strength data are shown in Figures 4-9.The estimates of PPDs for each of the input parameters required by the fluid, grain-shearing elastic, and poroelastic models were obtained by kernel density estimation based on the DREAM samples.

Fluid Model
Figure 4 displays the estimates of one-dimensional marginal PPDs for the fluid model parameters.Inspection of the figure reveals that the roughness spectral exponent γ 2 , roughness spectral strength ω 2 , density fluctuation spectral exponent γ 3 , and porosity β were well resolved by the inversion, as their PPDs are much narrower than their respective priors.This is consistent with the results of the sensitivity analysis.The density fluctuation spectral strength ω 3 , the ratio of compressibility to density fluctuations µ, the ratio of mass density of grains to water ρ r , and the ratio of bulk modulus of grains to water K r were somewhat resolved by the inversion, as their PPDs changed appreciably from their respective priors.This is not clearly seen in the results of the sensitivity analysis due to its limited resolution.The PPDs of the remaining parameters (mean grain diameter d, tortuosity α, permeability κ) only have a weak peak compared to their priors and provide an unreliable parameter inference.

Parameter Inversion
The parameter priors used for the parameter inversion were the same as those for the sensitivity analysis as shown in Tables 1 and 2. A total of five chains were used during the sampling process and deemed sufficient for parameter dimensions of all models.Under the premise of final convergence, the number of total generations varied between 120,000 and 400,000 depending on the model.The PPDs and Bayes factors were estimated by Equations ( 51), (52), and (61) using the last 60,000 DREAM samples.Estimates of the one-and two-dimensional marginal PPDs for the three models obtained via Bayesian parameter inversion applied to the backscattering strength data are shown in Figures 4-9.The estimates of PPDs for each of the input parameters required by the fluid, grain-shearing elastic, and poroelastic models were obtained by kernel density estimation based on the DREAM samples.

Fluid Model
Figure 4 displays the estimates of one-dimensional marginal PPDs for the fluid model parameters.Inspection of the figure reveals that the roughness spectral exponent  , roughness spectral strength  , density fluctuation spectral exponent  , and porosity  were well resolved by the inversion, as their PPDs are much narrower than their respective priors.This is consistent with the results of the sensitivity analysis.The density fluctuation spectral strength  , the ratio of compressibility to density fluctuations , the ratio of mass density of grains to water  , and the ratio of bulk modulus of grains to water  were somewhat resolved by the inversion, as their PPDs changed appreciably from their respective priors.This is not clearly seen in the results of the sensitivity analysis due to its limited resolution.The PPDs of the remaining parameters (mean grain diameter , tortuosity , permeability ) only have a weak peak compared to their priors and provide an unreliable parameter inference.The two-dimensional characteristics obtained by the inversion of the fluid model are presented in Figure 5, where the two-dimensional marginal PPDs are shown on the top right, and parameter correlations are shown on the bottom left.For intuitive and qualitative understanding, the two-dimensional PPDs have their own color bars, which are not given in the figure.The quantitative correlations are given by the correlation matrix map.The same display applies to the two-dimensional distribution for the next two models.It can be seen that parameter pairs with appreciable positive correlations are as follows: roughness spectral exponent γ 2 and roughness spectral strength ω 2 , roughness spectral exponent γ 2 and porosity β, roughness spectral strength ω 2 and porosity β, and the ratio of bulk modulus of grains to water K r and porosity β.Parameter pairs with appreciable negative correlations are the ratio of compressibility to density fluctuations µ and density fluctuations spectral exponent γ 3 .Note that parameters with obvious correlation are basically the parameters that are well resolved, as shown in Figure 4.When compared with the roughness spectral exponent γ 2 , the roughness spectral strength ω 2 with a relatively low Borgonovo index obtains satisfactory resolution by inversion due to the correlation, and this also appears in the next two inversion results.

Grain-Shearing Elastic Model
Estimates of the one-dimensional marginal PPDs for the grain-shearing elastic model parameters obtained via inversion applied to backscattering strength are shown in Figure 6.We can see by comparing PPDs and priors that the roughness spectral exponent  , roughness spectral strength  , density fluctuation spectral exponent  , and porosity  are resolved to a relatively high degree by the inversion.It can be noticed that the resolution of porosity  for the grain-shearing elastic model is weaker than that for the fluid model.This could be due to the difference in density ratios used in

Grain-Shearing Elastic Model
Estimates of the one-dimensional marginal PPDs for the grain-shearing elastic model parameters obtained via inversion applied to backscattering strength are shown in Figure 6.We can see by comparing PPDs and priors that the roughness spectral exponent γ 2 , roughness spectral strength ω 2 , density fluctuation spectral exponent γ 3 , and porosity β are resolved to a relatively high degree by the inversion.It can be noticed that the resolution of porosity β for the grain-shearing elastic model is weaker than that for the fluid model.This could be due to the difference in density ratios used in the two models, as shown in Equations ( 33) and ( 17), while sediment mass density is most closely related to porosity.The density fluctuation spectral strength ω 3 , the ratio of compressibility to density fluctuation µ, the ratio of mass density of grains to water ρ r , the ratio of bulk modulus of grains to water K r , material exponent n, and compressional viscoelastic relaxation time τ 1 are somewhat resolved by the inversion.The PPDs of the compressional rigidity coefficient γ 1 and shear rigidity coefficient γ t did not appreciably change from their priors.It is not difficult to find that the results in Figure 6 are consistent with most of the previous results of sensitivity analysis in Figure 1b.rigidity coefficient  did not appreciably change from their priors.It is not difficult to find that the results in Figure 6 are consistent with most of the previous results of sensitivity analysis in Figure 1b.Included in Figure 7 are the estimated two-dimensional marginal PPDs and parameter correlations obtained by the inversion based on the grain-shearing elastic model.From Figure 7 we can see that the parameter correlations of the grain-shearing elastic model are similar to those for the fluid model.In addition, the unique parameter pair of grain-shearing elastic models with appreciable positive correlation is the ratio of bulk modulus of grains to water  and material exponent .The parameter pair with appreciable negative correlation is compressional rigidity coefficient  and material exponent .Included in Figure 7 are the estimated two-dimensional marginal PPDs and parameter correlations obtained by the inversion based on the grain-shearing elastic model.From Figure 7 we can see that the parameter correlations of the grain-shearing elastic model are similar to those for the fluid model.In addition, the unique parameter pair of grain-shearing elastic models with appreciable positive correlation is the ratio of bulk modulus of grains to water K r and material exponent n.The parameter pair with appreciable negative correlation is compressional rigidity coefficient γ 1 and material exponent n.
Included in Figure 7 are the estimated two-dimensional marginal PPDs and parameter correlations obtained by the inversion based on the grain-shearing elastic model.From Figure 7 we can see that the parameter correlations of the grain-shearing elastic model are similar to those for the fluid model.In addition, the unique parameter pair of grain-shearing elastic models with appreciable positive correlation is the ratio of bulk modulus of grains to water  and material exponent .The parameter pair with appreciable negative correlation is compressional rigidity coefficient  and material exponent .

Poroelastic Model
The estimated one-dimensional marginal PPDs for the parameters of the poroelastic model are shown in Figure 8. Comparing the marginal PPDs to their respective priors indicates that the roughness spectral exponent γ 2 , roughness spectral strength ω 2 , density fluctuation spectral exponent γ 3 , ratio of compressibility to density fluctuation µ, porosity β, ratio of mass density of grains to water ρ r , ratio of bulk modulus of grains to water K r , and cementation exponent m are well resolved by the inversion.The resolution of porosity β of the poroelastic model is similar to that of the grain-shearing elastic model, and this may be due to the same form of density ratio.The density fluctuation spectral strength ω 3 and high-frequency asymptotic increase K y are somewhat resolved by the inversion, and the others (mean grain diameter d, pore shape factor a B , Poisson's ratio of grains ν, low-frequency asymptotic frame bulk modulus K bo , bulk relaxation frequency f k ) have unreliable resolutions.The resolutions of the two ratios of the poroelastic model, ρ r and K r , are slightly better than those of the previous two models.The possible reason is that the poroelastic model treats both porosity and elasticity and has inherent advantages for porous mediums such as sandy sediments, for which it is possible that the fluid and granular phases will vibrate differently in response to acoustic excitation [3].From a mathematical perspective, the two parameters can be more coupled with the scattering strength, as shown by the equations in Appendix B. The resolutions of interface roughness parameters (roughness spectral exponent γ 2 , roughness spectral strength ω 2 ) and volume inhomogeneity parameters (density fluctuation spectral exponent γ 3 , density fluctuation spectral strength ω 3 , ratio of compressibility to density fluctuations µ) of the three models have a certain degree of difference, while the common parameters of the poroelastic model have slightly better resolution in general.Compared with the previous sensitivity analysis results in Figure 1c, the results here show better resolution and have more information about model parameters.Despite the total number of model parameters and the number of parameters that are somewhat resolved, the inversion based on the poroelastic model provides four more well-resolved parameters than that based on the other two models.The resolutions of common parameters of the poroelastic model are also slightly better than that of the other two models.From the perspective of the PPDs of the model parameters, more reliable information of sediment properties is provided by the inversion based on the poroelastic model, which can be considered superior to the other two models.Figure 9 shows the estimates of two-dimensional marginal PPDs and parameter correlations obtained by the inversion based on the poroelastic model.In addition to the parameter correlations similar to the fluid model, the parameter pair, which is the ratio of bulk modulus of grains to water  and cementation exponent , has an appreciable positive correlation.It should be viewed with caution that these correlations are only shown by the samples, and not all of them have practical physical meaning.The same attention needs to be paid to the two-dimensional marginal PPDs of the previous two models.Indeed, it is difficult to provide a physical reason for the parameter correlations shown by the inversion results, as part of them cannot be measured by any means.The authors believe that the universality of these correlations needs to be further verified and analyzed by modelbased inversions using additional measurements for different types of sandy sediments.Taking available sediment properties measured in situ as intermediate quantities of parameter pairs with correlations may help explain the correlations.Figure 9 shows the estimates of two-dimensional marginal PPDs and parameter correlations obtained by the inversion based on the poroelastic model.In addition to the parameter correlations similar to the fluid model, the parameter pair, which is the ratio of bulk modulus of grains to water K r and cementation exponent m, has an appreciable positive correlation.It should be viewed with caution that these correlations are only shown by the samples, and not all of them have practical physical meaning.The same attention needs to be paid to the two-dimensional marginal PPDs of the previous two models.Indeed, it is difficult to provide a physical reason for the parameter correlations shown by the inversion results, as part of them cannot be measured by any means.The authors believe that the universality of these correlations needs to be further verified and analyzed by model-based inversions using additional measurements for different types of sandy sediments.Taking available sediment properties measured in situ as intermediate quantities of parameter pairs with correlations may help explain the correlations.

Model Comparison and Selection
The model comparison and selection results obtained by the inversion using backscattering strength data are shown in Figures 10-12.Figure 10 gives a comparison of the uncertainties of the ratios of compressional and shear sound speed in sediment to compressional sound speed in water.The uncertainties of attenuation of the two types of sound waves are also shown in Figure 10.The sets of estimated sound speed ratios and attenuation at different frequencies are calculated by the samples with the three geoacoustic models, EDFM, VGS(λ), and CREB.The PPDs of sound speed ratio and attenuation are also obtained by kernel density estimation based on the estimated sets.Because of the finite number of three frequencies, the frequency dependences for sound speed and attenuation are not obvious in these results.Comparing the PPDs of compressional wave indicates that the predicted uncertainties of CREB are slightly more concentrated than those of EDFM and their maximum posterior estimates are very close to each other, as shown in the left side of Figure 10.The PPDs of the compressional wave of EDFM and CREB have a slight sound speed ratio difference and an appreciable attenuation difference compared to VGS(λ).This suggests that the scattering strengths of the fluid and poroelastic models are more sensitive to attenuation than the grain-shearing elastic model.Besides, the scattering strengths of the three models are much more sensitive to compressional sound speed than attenuation as the PPDs of compressional speed ratio are all highly concentrated.

Model Comparison and Selection
The model comparison and selection results obtained by the inversion using backscattering strength data are shown in Figures 10-12.Figure 10 gives a comparison of the uncertainties of the ratios of compressional and shear sound speed in sediment to compressional sound speed in water.The uncertainties of attenuation of the two types of sound waves are also shown in Figure 10.The sets of estimated sound speed ratios and attenuation at different frequencies are calculated by the samples with the three geoacoustic models, EDFM, VGS(λ), and CREB.The PPDs of sound speed ratio and attenuation are also obtained by kernel density estimation based on the estimated sets.Because of the finite number of three frequencies, the frequency dependences for sound speed and attenuation are not obvious in these results.Comparing the PPDs of compressional wave indicates that the predicted uncertainties of CREB are slightly more concentrated than those of EDFM and their maximum posterior estimates are very close to each other, as shown in the left side of Figure 10.The PPDs of the compressional wave of EDFM and CREB have a slight sound speed ratio difference and an appreciable attenuation difference compared to VGS(λ).This suggests that the scattering strengths of the fluid and poroelastic models are more sensitive to attenuation than the grain-shearing elastic model.Besides, the scattering strengths of the three models are much more sensitive to compressional sound speed than attenuation as the PPDs of compressional speed ratio are all highly concentrated.From the right side of Figure 10, it can be seen that shear sound speed predicted by CREB is larger than that predicted by VGS( λ ) and the uncertainties of attenuation are roughly similar.Compared with compressional sound speed, the shear sound speed in sandy sediments is much lower.This suggests that the addition of shear waves is only a minor modification to model the scattering strength [3] and the uncertainty of shear sound speed is larger than that of compressional sound speed.Since slow compressional waves supported in a poroelastic medium are given theoretically and have never been observed in natural marine sediment, the uncertainties of sound speed and attenuation for slow compressional waves are not given here.
Figure 11 depicts a comparison of the uncertainties of backscattering strength predicted by the samples.The DREAM samples were substituted into each model to calculate the estimated backscattering strength set at different frequencies and angles.For a selected model at a selected frequency, the estimated backscattering strength set at each angle was statistically transformed into a probability distribution, which is represented by a vertical color strip.All the vertical color strips for the whole angle range were spliced together along the x-axis to form the subgraphs in Figure 11.Experimental measurements are also shown in the figure.In general, all three models are well matched with the experimental measurements.Appreciable differences of uncertainties of the estimated backscattering strength for the three models can be seen in Figure 11, and the fluid model has the most concentrated posterior predictive distributions.The results of the fluid model are only slightly different from those of the poroelastic model near the critical angle.Compared with the fluid and poroelastic models, the results for the grain-shearing elastic model have the greatest uncertainties and smallest concentration.The fluid model seemingly does the best job of matching the frequency and angle dependence of the backscattering strength data for sandy sediments.Comparing Figures 10 and 11, it can also be seen that differences of uncertainties of the estimated backscattering strength correspond to the differences of uncertainties of the estimated sound speed and attenuation for the From the right side of Figure 10, it can be seen that shear sound speed predicted by CREB is larger than that predicted by VGS(λ) and the uncertainties of attenuation are roughly similar.Compared with compressional sound speed, the shear sound speed in sandy sediments is much lower.This suggests that the addition of shear waves is only a minor modification to model the scattering strength [3] and the uncertainty of shear sound speed is larger than that of compressional sound speed.Since slow compressional waves supported in a poroelastic medium are given theoretically and have never been observed in natural marine sediment, the uncertainties of sound speed and attenuation for slow compressional waves are not given here.
Figure 11 depicts a comparison of the uncertainties of backscattering strength predicted by the samples.The DREAM samples were substituted into each model to calculate the estimated backscattering strength set at different frequencies and angles.For a selected model at a selected frequency, the estimated backscattering strength set at each angle was statistically transformed into a probability distribution, which is represented by a vertical color strip.All the vertical color strips for the whole angle range were spliced together along the x-axis to form the subgraphs in Figure 11.Experimental measurements are also shown in the figure.In general, all three models are well matched with the experimental measurements.Appreciable differences of uncertainties of the estimated backscattering strength for the three models can be seen in Figure 11, and the fluid model has the most concentrated posterior predictive distributions.The results of the fluid model are only slightly different from those of the poroelastic model near the critical angle.Compared with the fluid and poroelastic models, the results for the grain-shearing elastic model have the greatest uncertainties and smallest concentration.The fluid model seemingly does the best job of matching the frequency and angle dependence of the backscattering strength data for sandy sediments.Comparing Figures 10  and 11, it can also be seen that differences of uncertainties of the estimated backscattering strength correspond to the differences of uncertainties of the estimated sound speed and attenuation for the three models.The fluid and poroelastic models are similar, and they both have appreciable differences compared with the grain-shearing elastic model.[10,12,19].The poroelastic model, which takes sediment composed of grains and water as a two-phase system, is essentially the most reasonable approximation for sandy sediments.In this sense, the poroelastic and fluid models are both superior to the grain-shearing elastic model.Although the Bayes factors show slight differences between the fluid and poroelastic models, we still consider the poroelastic model to be the best because the inversion based on the poroelastic model provides the best-resolved parameters, as discussed in Section 4.2.3.[10,12,19].The poroelastic model, which takes sediment composed of grains and water as a two-phase system, is essentially the most reasonable approximation for sandy sediments.In this sense, the poroelastic and fluid models are both superior to the grain-shearing elastic model.Although the Bayes factors show slight differences between the fluid and poroelastic models, we still consider the poroelastic model to be the best because the inversion based on the poroelastic model provides the best-resolved parameters, as discussed in Section 4.2.3.

Conclusions
Based on the combination of geoacoustic models and corresponding scattering models, Bayesian inversion and model selection techniques are applied to compare three sediment acoustic models: the fluid model with EDFM, the grain-shearing elastic model with VGS(λ), and the poroelastic model with CREB.The combination allows us to use backscattering strength measured remotely to carry out the parameter inversions.Parameter uncertainties are given by the one-and two-dimensional PPDs estimated by the DREAM samples.
Two studies were carried out using the inversion results.First, the resolution and correlation of parameters for the three models were compared based on the estimates of PPDs.The poroelastic model provides more well-resolved parameters than the other two models.In general, the resolutions of common parameters of the poroelastic model are also slightly better than those of the other two models.In other words, the inversion based on the poroelastic model can provide more highly resolved information about sediments properties, some of which are difficult to measure directly.
The second study involved model comparison and selection using model predictions and Bayes factors.Based on the DREAM samples, the predictions using different models and acoustic frequencies were calculated for qualitative comparative analysis with the backscattering strength measurements, and the Bayes factors of the three models are calculated for quantitative comparison.From the analysis, it appears that all three models can produce predictions that agree with the measurements, though the fluid and poroelastic models appear to outperform the grain-shearing elastic model.It should be viewed with caution that the model comparison results are obtained by the inversion based only on the limited measurements of sandy sediments.Indeed, the effectiveness of the poroelastic model for sandy sediments has been authenticated and is noticeable.
The results in this paper indicate that the best model for sandy sediments is the poroelastic model.However, it has been suggested that the elastic model presumably is better for softer sediments than the poroelastic model.Besides, the interface roughness scattering model that accommodates shear can provide an appreciable correction to the fluid model for consolidated sediments such as rock seafloor, which has a greater shear wave speed than the sound speed in water.Much attention has been paid to model sound propagation in sandy sediments, and no single model can be applied to all types of sediment.Despite its importance, this topic has largely been ignored in the literature to date.We have performed a comparative study of evidence estimation for the remote sensing of sediment properties in this paper.This was an attempt to implement a combination of several different models for geoacoustic inversion, and the model that best matches the

Conclusions
Based on the combination of geoacoustic models and corresponding scattering models, Bayesian inversion and model selection techniques are applied to compare three sediment acoustic models: the fluid model with EDFM, the grain-shearing elastic model with VGS(λ), and the poroelastic model with CREB.The combination allows us to use backscattering strength measured remotely to carry out the parameter inversions.Parameter uncertainties are given by the one-and two-dimensional PPDs estimated by the DREAM samples.
Two studies were carried out using the inversion results.First, the resolution and correlation of parameters for the three models were compared based on the estimates of PPDs.The poroelastic model provides more well-resolved parameters than the other two models.In general, the resolutions of common parameters of the poroelastic model are also slightly better than those of the other two models.In other words, the inversion based on the poroelastic model can provide more highly resolved information about sediments properties, some of which are difficult to measure directly.
The second study involved model comparison and selection using model predictions and Bayes factors.Based on the DREAM samples, the predictions using different models and acoustic frequencies were calculated for qualitative comparative analysis with the backscattering strength measurements, and the Bayes factors of the three models are calculated for quantitative comparison.From the analysis, it appears that all three models can produce predictions that agree with the measurements, though the fluid and poroelastic models appear to outperform the grain-shearing elastic model.It should be viewed with caution that the model comparison results are obtained by the inversion based only on the limited measurements of sandy sediments.Indeed, the effectiveness of the poroelastic model for sandy sediments has been authenticated and is noticeable.
The results in this paper indicate that the best model for sandy sediments is the poroelastic model.However, it has been suggested that the elastic model presumably is better for softer sediments than the poroelastic model.Besides, the interface roughness scattering model that accommodates shear can provide an appreciable correction to the fluid model for consolidated sediments such as rock seafloor, which has a greater shear wave speed than the sound speed in water.Much attention has been paid to model sound propagation in sandy sediments, and no single model can be applied to all types of sediment.Despite its importance, this topic has largely been ignored in the literature to date.We have performed a comparative study of evidence estimation for the remote sensing of sediment properties in this paper.This was an attempt to implement a combination of several different models for geoacoustic inversion, and the model that best matches the measurements was obtained.The attempt provides a relatively feasible remote sensing scheme for various types of sediments under unknown conditions.However, this scheme should be regarded as tentative and requires more measurements of different types of sediment for validation.where the vertical line in Equation(A20) separates it into the 5 × 5 matrices,  ( ) () , and the column matrices,  ( ) ().The elements of  ( ) () and  ( ) () are: where the vertical line in Equation (A20) separates it into the 5 × 5 matrices, P (n) (K), and the column matrices, Q (n) (K).The elements of P (n) (K) and Q (n) (K) are:

Figure 3
Figure3shows measured backscattering strength as a function of grazing angle at different frequencies.It can be seen that the backscattering strength has only a weak frequency dependence.There is a slight drop of backscattering strength above the critical angle of 30° which is consistent with the measurements given by Williams[10].Although no field experiments were carried out, the authors believe that laboratory measurements are more manageable and provide adequate fitness between model and data.Of course, more subsequent field data for model selection analysis are required due to the limitations of laboratory measurements.This work features the use of a remote sensing scheme with multiple models and lays a foundation for future practical field applications.

Figure 3
Figure3shows measured backscattering strength as a function of grazing angle at different frequencies.It can be seen that the backscattering strength has only a weak frequency dependence.There is a slight drop of backscattering strength above the critical angle of 30 • which is consistent with the measurements given by Williams[10].Although no field experiments were carried out, the authors believe that laboratory measurements are more manageable and provide adequate fitness between model and data.Of course, more subsequent field data for model selection analysis are required due to the limitations of laboratory measurements.This work features the use of a remote sensing scheme with multiple models and lays a foundation for future practical field applications.

Figure 3 .
Figure 3. Measured backscattering strength as a function of grazing angle of different frequencies.

30 Figure 3 .
Figure 3. Measured backscattering strength as a function of grazing angle of different frequencies.

Figure 4 .Figure 4 .
Figure 4. Estimates of one-dimensional marginal posterior probability distributions (PPDs) for fluid model parameters.The two-dimensional characteristics obtained by the inversion of the fluid model are presented in Figure 5, where the two-dimensional marginal PPDs are shown on the top right, and parameter correlations are shown on the bottom left.For intuitive and qualitative understanding, the twodimensional PPDs have their own color bars, which are not given in the figure.The quantitative Remote Sens. 2018, 10, x FOR PEER REVIEW 17 of 30 appreciable negative correlations are the ratio of compressibility to density fluctuations  and density fluctuations spectral exponent  .Note that parameters with obvious correlation are basically the parameters that are well resolved, as shown in Figure 4.When compared with the roughness spectral exponent  , the roughness spectral strength  with a relatively low Borgonovo index obtains satisfactory resolution by inversion due to the correlation, and this also appears in the next two inversion results.

Figure 5 .
Figure 5. Estimates of two-dimensional marginal PPDs for fluid model parameters (top right) and correlation matrix map (bottom left).

Figure 5 .
Figure 5. Estimates of two-dimensional marginal PPDs for fluid model parameters (top right) and correlation matrix map (bottom left).

Figure 6 .
Figure 6.Estimates of one-dimensional marginal PPDs for grain-shearing elastic model parameters.

Figure 6 .
Figure 6.Estimates of one-dimensional marginal PPDs for grain-shearing elastic model parameters.

Figure 7 .
Figure 7. Estimates of two-dimensional marginal PPDs for grain-shearing elastic model parameters (top right) and correlation matrix map (bottom left).

Figure 8 .
Figure 8. Estimates of one-dimensional marginal PPDs for poroelastic model parameters.

Figure 9 .
Figure 9. Estimates of two-dimensional marginal PPDs for poroelastic model parameters (top right) and correlation matrix map (bottom left).

Figure 9 .
Figure 9. Estimates of two-dimensional marginal PPDs for poroelastic model parameters (top right) and correlation matrix map (bottom left).

Figure 10 .
Figure 10.Comparison of uncertainties of ratios of compressional and shear sound speed in sediment to compressional sound speed in water and corresponding attenuation.

Figure 10 .
Figure 10.Comparison of uncertainties of ratios of compressional and shear sound speed in sediment to compressional sound speed in water and corresponding attenuation.
Remote Sens. 2018, 10, x FOR PEER REVIEW 23 of 30 three models.The fluid and poroelastic models are similar, and they both have appreciable differences compared with the grain-shearing elastic model.

Figure 11 .
Figure 11.Comparison of uncertainties of backscattering strength predicted by samples based on different models and acoustic frequencies.Red asterisks represent corresponding experimental measurements, as shown in Figure 3.

Figure 11 .
Figure 11.Comparison of uncertainties of backscattering strength predicted by samples based on different models and acoustic frequencies.Red asterisks represent corresponding experimental measurements, as shown in Figure 3.The comparison of Bayes factors of different model pairs obtained by the inversion based on the backscattering strength data is summarized in Figure 12.Referring to the factor scale given in Table 3, the evidence in favor of the fluid model over the grain-shearing elastic model can be considered substantial.The Bayesian evidence for the fluid and poroelastic models is similar, with the fluid model slightly favored.This similarity and the similarity of compressional sound speed and attenuation, as shown in Figure 10, can be understood since the EDFM used in the fluid model is essentially an approximation of the full Biot-Stoll model and the corrected and reparametrized extended Biot-Stoll model is used in the poroelastic model.It has been proven that backscattering strength, compressional sound speed, and attenuation predicted by the fluid model based on the EDFM are very close to the predictions of full Biot-Stoll theory for sandy sediment[10,12,19].The poroelastic model, which takes sediment composed of grains and water as a two-phase system, is essentially the most reasonable approximation for sandy sediments.In this sense, the poroelastic and fluid models are both superior to the grain-shearing elastic model.Although the Bayes factors show slight differences between the fluid and poroelastic models, we still consider the poroelastic model to be the best because the inversion based on the poroelastic model provides the best-resolved parameters, as discussed in Section 4.2.3.

Figure 12 .
Figure 12.Comparison of Bayes factors calculated for each pair of models through Bayesian inference applied to backscattering strength data.Number superimposed over each matrix element is the base 10 logarithm of the Bayes factor for a given pair of models.

Figure 12 .
Figure 12.Comparison of Bayes factors calculated for each pair of models through Bayesian inference applied to backscattering strength data.Number superimposed over each matrix element is the base 10 logarithm of the Bayes factor for a given pair of models.

Table 1 .
Parameter bounds used to construct uniform priors for models' input parameters.

Table 2 .
Known parameters for water.

Table 3 .
Bayes factor scale and corresponding interpretation.