A Passive Source Location Method in a Shallow Water Waveguide with a Single Sensor Based on Bayesian Theory

Bayesian methodology is a good way to infer unknown parameters in a marine environment. A passive source location method in a shallow water waveguide with a single sensor based on Bayesian theory is presented in this paper. The input of a Bayesian inversion algorithm is received different normal mode impulse signals, which are separated and extracted with a warping transformation from received broadband impulse signals. The source range, depth, and other seabed parameters were estimated without prior knowledge of the seabed information. Different normal mode impulse acoustic signals travelling at different group speeds arrived at the sensor at different times because of the dispersion characteristics of the shallow water waveguide. The time delay of different modes can be used for the passive source location. However, normal mode group speeds are greatly affected by the environmental parameters. The performance of the passive location becomes negative when parameters mismatch. In this paper, the source location was transformed to the inversion of the source location and environmental parameters, which can be estimated accurately based on the multi-dimensional posterior probability density (PPD). This method is less limited by environmental factors, and the accuracy of inversion results can be analyzed according to the PPD of inversion parameters, which has higher reliability and a wider application scope. The effectiveness and robustness of the algorithm were quantified in terms of the root mean squared error (RMSE) at a variety of signal-to-noise ratios (SNRs) in 50 simulation sets. The RMSE values decreased with the SNR. The validity and accuracy of the method were proved by the results of simulation and experiment data.


Introduction
The passive location of underwater sources is the basis of underwater target detection and is an active field in underwater acoustic research. Based on the number of sensors, passive source location technology in a shallow water waveguide can be divided into two categories: methods using a sensor array and those with a single sensor.
For methods with a sensor array (vertical array or horizontal array), the water column is covered by sensors, and the marine environment information can be measured after processing the received signal data. The source will be estimated accurately using the marine environment information. However, in practical applications, a sensor array carries a high computational cost; additionally, called the genetic algorithm (GA) is used to search the optimum solution [23], and the source range and depth are then inverted. The method is applied to measured data collected in a shallow water experiment in 2014, and the results compare well with global positioning system (GPS) measurements taken during the experiment.
The remainder of this paper is organized as follows: Section 2 describes the separation process for the received signal and the extraction of normal mode signals. Section 3 presents the location theory and algorithms, while Section 4 applies the location method to the simulated data. Section 5 presents and discusses the location results of the experimental data. Finally, Section 6 provides concluding remarks.

Modal Propagation Theory
In normal mode theory, the received pressure field signal is the sum of several normal mode signals. Assuming a broadband source emitting at depth z s in a range-independent shallow water waveguide with a half-infinite liquid seabed, the pressure field Y(f,r) received at depth z r after propagation over a range r can be written as [24] where S(f ) is the source spectrum, N is the number of normal modes, ξ n (f ) is the horizontal wavenumber, and U n is the modal depth function, whose amplitude varies with depth, β n (f ) is the attenuation coefficient, which is a very small value and varies with frequency, and ρ(z s ) is the density at the source depth z s , which is the density of water.

Warping
Transformation Y(f,r) consists of several normal mode signals of different numbers, each number signal contains a significant amount of environmental information, so that a single normal mode signal can also be used for passive source location and parameter inversion after analysis and processing.
A signal processing method called warping [25] is used to separate and extract the normal mode signals from the received pressure field signal. The warping transformation is a model-based transformation designed to linearize the signal phase. Each mode will become a single-frequency signal with its invariant frequency following the warping transformation.
For a given signal y(t), which is the received signal in the time domain, the warping transformation can be written as [21] where Wy(t) is the warped signal, and h(t) is the warping operator while ∂t is the derivative of w(t). Warping transformation is reversible; the warped signal can be unwarped using h −1 (t), and |∂h(t)/∂t| only provides the energy conservation between the original and warped signals.
The warped operator can be written as and the unwarped operator is where t r = r c 0 , and c 0 is the water sound speed. Warping transformation is robust and can be applied to most low frequency shallow water scenarios without detailed knowledge of the environment or precise propagation range [26]. Therefore, t r can be determined empirically without knowing r or c 0 . The received signal can be presented in the TF domain using short-term Fourier transformation (STFT) [27]. Warping transformation is then used in the TF domain, the received signal becomes the sum of several number linear single frequency modes, and the correspondent frequency of each mode is the eigen-frequency. The single normal mode signal is extracted by a frequency filter. The signal in the frequency domain can be transformed from the time domain using Fourier transformation. Therefore, when warping transformation is used, the n number normal mode pressure field signal at frequency f can be written as The procedure for the separation and extraction of the normal mode signals is described in Figure 1.  [27]. Warping transformation is then used in the TF domain, the received signal becomes the sum of several number linear single frequency modes, and the correspondent frequency of each mode is the eigen-frequency. The single normal mode signal is extracted by a frequency filter. The signal in the frequency domain can be transformed from the time domain using Fourier transformation. Therefore, when warping transformation is used, the n number normal mode pressure field signal at frequency f can be written as The procedure for the separation and extraction of the normal mode signals is described in Figure 1.

Bayesian Inversion Theory
In a Bayesian inversion, the multi-dimensional posterior probability density (PPD) is usually interpreted in terms of model-parameter estimates and uncertainties. In a Bayesian approach, let m be a vector of M free parameters representing a realization, and let d represent N measured data which constrain the model. These quantities are considered random variables that are related via Bayes' rule [28]: The P(m|d) is the PPD, P(m) is the prior distribution, P(d) is the probability density of the measured data, which is independent of m, and P(d|m) represents the conditional probability density for d, which is interpreted as the likelihood L(m) for the measured data. Thus, Equation (6) can be written as where The likelihood function depends on the statistical distribution of the data expression and errors (measurement error and theoretical error), and E(m) is the data misfit function. In cases where the error distribution is not known independently, a good strategy is to choose the Gaussian distribution

Bayesian Inversion Theory
In a Bayesian inversion, the multi-dimensional posterior probability density (PPD) is usually interpreted in terms of model-parameter estimates and uncertainties. In a Bayesian approach, let m be a vector of M free parameters representing a realization, and let d represent N measured data which constrain the model. These quantities are considered random variables that are related via Bayes' rule [28]: P(m |d ) = P(d |m )P(m)/P(d).
The P(m|d) is the PPD, P(m) is the prior distribution, P(d) is the probability density of the measured data, which is independent of m, and P(d|m) represents the conditional probability density for d, which is interpreted as the likelihood L(m) for the measured data. Thus, Equation (6) can be written as where The likelihood function depends on the statistical distribution of the data expression and errors (measurement error and theoretical error), and E(m) is the data misfit function. In cases where the error distribution is not known independently, a good strategy is to choose the Gaussian distribution and estimate the statistical parameters from the data. A generalized misfit combining data and prior can be defined as where ϕ(m) is the cost function. PPD is written as The integration domain spans an M-dimensional parameter space. In Bayesian inversion, the PPD of m is interpreted as the inversion results. In this paper, the maximum a posteriori (MAP) model is used, and the expression iŝ Based on Equation (9), Equation (11) can be written aŝ Additionally, the mean model, the posterior model covariance matrix, and the one-and two-dimensional marginal probability densities are defined respectively as where (m − E(m )) T is the transpose of (m − E(m )).

Cost Function
To estimate the marine environment parameters using Bayesian inversion methodology, a sufficient cost function is necessary.
The traditional Bayesian inversion method is usually carried out by a sensor array (horizontal array or vertical array). Not only can the accuracy of measured data be affected by the array form, the array also carries a significant computational cost when used. In this paper, the different number normal mode signals in frequency are the input of Bayesian methodology, and estimation of the source range and depth is carried out by matching different amounts of normal mode signals in the frequency domain.
Consider the measured data P m , which is a matrix of N × N f , where N is the number of modes, and N f is the number of frequency points. The element of P m at line n and column n f is therefore where ∑ n=1 n = N and ∑ n f =1 n f = N f . Assuming the data errors are Gaussian-distributed random variables and covariance matrix C n f , the likelihood function is where P m n f represents the measured data at the n f -th frequency point of P c , P c n f represents the modeled data at the n f -th frequency point of P c , and P c can be written as P c (m) = H(m)S (19) where H(m) is the channel transition function, and S is the source spectrogram. In many cases, the error statistics are unknown; the covariance matrices C n f should be estimated from data. Consider first the common approximation of the independent, identical errors and diagonal covariance matrices C n f = υ 2 n f I, where υ 2 n f is the variance for the n f -th frequency point and I is the identity matrix. The likelihood function is where S nf is the source spectrum for the n f -th frequency point. The likelihood function L(m) can be expressed as the experiential index of the cost function ϕ(m), so the likelihood function is The cost function is The υ n f and S can be estimated by maximizing the likelihood [29], setting (24) and Equation (24) can be written as where H T n f is the conjugate transpose matrix of H n f , and the estimation of the source spectrum is The S n f is substituted for S n f in Equation (20), and the cost function can be written as According to ∂ϕ ∂υ n f = 0, the estimation of υ 2 n f is Therefore, the cost function is When the constant term is ignored, the cost function is The source range, depth, and other seabed parameters can be inverted by minimizing ϕ(m). The GA is used to search for the optimum solution, when the optimal value of search does not change, and converges to a fixed value, which is considered the optimal value. The mutation probability, selection probability, crossover probability, generation number, and initial population size are 0.05, 0.5, 0.8, 5000, and 64, respectively. In addition, 20 sets are computed in parallel to ensure the parameters converge to the global optimum.

Simulation Example
The simulation was performed in a shallow water waveguide with a half-infinite liquid seabed. The depth was 25 m, the sound speed in water was an isovelocity with c 0 = 1500 m/s, a broadband source was emitted at depth z s = 20 m with frequency band 200~300 Hz, the source is a linear frequency modulated impulse signal, the SNR was 20 dB, and the signal was received at depth z r = 23 m after propagation range r = 7700 m. The seabed sound speed was c b = 1650 m/s, and seabed density was

The Extraction of Normal Mode
The received signal in the time domain and the extracted modes are shown in Figure 2. The received signal contains several normal modes from Figure 2b, Figure 2c shows that the normal mode can be separated after warping transformation, and the phases of the mode signals are transformed from non-linear frequency modulation to linear. Each mode becomes a single-frequency signal with its invariant frequency after warping transformation, and the spectrum characteristics are shown in Figure 2d. Warping transformation is reversible. The mode signals are extracted by a frequency filter and unwarping transformation. The extracted signals in time and TF domains of the first four modes are shown in Figure 3. Additionally, a comparison between the original received signal in the time domain and the signal recovered from the first four warped mode signals was made, and the result illustrates that the recovery signal is consistent with the original signal, but a small amount of noise was ignored ( Figure 4). Sensors 2019, 19, x FOR PEER REVIEW 8 of 18  (c) (d)  The input signals of the Bayesian inversion theory are the extracted normal mode signals in the frequency domain. The signal in the frequency domain can be transformed from the time domain using Fourier transformation. When the input signals are obtained, the source depth and range can be estimated based on Section 3.

The Analysis of the Inversion Parameter Sensitivity
To illustrate the validity of the cost function in Section 3.2, the inversion parameter sensitivity to the cost function is analyzed. When the parameter sensitivity is analyzed, the parameter to be analyzed is constantly changing in a certain range, while the other three parameters remain unchanged. When the analyzed parameter is near the true value, the cost function obtains the minimum value. Four parameters of source range, depth, seabed sound speed, and seabed density are analyzed. The results are as seen in Figure 5; parameters are sensitive to the cost function. The cost function is minimized at the true values, and the analysis curve of parameters varies sharply near the true value, so all inversion parameters are sensitive. The cost function is valid with respect to the inverse parameters.    The input signals of the Bayesian inversion theory are the extracted normal mode signals in the frequency domain. The signal in the frequency domain can be transformed from the time domain using Fourier transformation. When the input signals are obtained, the source depth and range can be estimated based on Section 3.

The Analysis of the Inversion Parameter Sensitivity
To illustrate the validity of the cost function in Section 3.2, the inversion parameter sensitivity to the cost function is analyzed. When the parameter sensitivity is analyzed, the parameter to be analyzed is constantly changing in a certain range, while the other three parameters remain unchanged. When the analyzed parameter is near the true value, the cost function obtains the minimum value. Four parameters of source range, depth, seabed sound speed, and seabed density are analyzed. The results are as seen in Figure 5; parameters are sensitive to the cost function. The cost function is minimized at the true values, and the analysis curve of parameters varies sharply near the true value, so all inversion parameters are sensitive. The cost function is valid with respect to the inverse parameters.  The input signals of the Bayesian inversion theory are the extracted normal mode signals in the frequency domain. The signal in the frequency domain can be transformed from the time domain using Fourier transformation. When the input signals are obtained, the source depth and range can be estimated based on Section 3.

The Analysis of the Inversion Parameter Sensitivity
To illustrate the validity of the cost function in Section 3.2, the inversion parameter sensitivity to the cost function is analyzed. When the parameter sensitivity is analyzed, the parameter to be analyzed is constantly changing in a certain range, while the other three parameters remain unchanged. When the analyzed parameter is near the true value, the cost function obtains the minimum value. Four parameters of source range, depth, seabed sound speed, and seabed density are analyzed. The results are as seen in Figure 5; parameters are sensitive to the cost function. The cost function is minimized at the true values, and the analysis curve of parameters varies sharply near the true value, so all inversion parameters are sensitive. The cost function is valid with respect to the inverse parameters.

The Inversion Results
The input signals of the inverse scheme are given in Section 4.1, and the replica (model signal) is computed by KRAKEN, an acoustic computation program [30]. All parameters were searched over relatively wide intervals based on Equation (30) by the GA, and the corresponding search bounds are given in Table 1. When the parameters of the search converge to the optimal value, the optimal value is placed into Equations (10) and (16), and the PPD and the marginal probability densities can be estimated. The inverse value of the parameter corresponds to the maximum marginal probability. Figure 6 presents the marginal probability densities for each individual parameter. 18 18.

The Inversion Results
The input signals of the inverse scheme are given in Section 4.1, and the replica (model signal) is computed by KRAKEN, an acoustic computation program [30]. All parameters were searched over relatively wide intervals based on Equation (30) by the GA, and the corresponding search bounds are given in Table 1. When the parameters of the search converge to the optimal value, the optimal value is placed into Equations (10) and (16), and the PPD and the marginal probability densities can be estimated. The inverse value of the parameter corresponds to the maximum marginal probability. Figure 6 presents the marginal probability densities for each individual parameter. To study the inter-relationship of these parameters, Figure 7 shows the correlation coefficient matrix of arbitrary two inversion parameters calculated by correlating the marginal probability densities of two parameters and normalizing them. Range r and source depth z s correlate with the seabed sound speed c b , which causes multiple peaks in PPD for c b . When the inversion parameters have multiple peaks in marginal probability densities, many sets computed in parallel are necessary to ensure that the parameters converge to the global optimum. Other parameters do not indicate a strong correlation. To study the inter-relationship of these parameters, Figure 7 shows the correlation coefficient matrix of arbitrary two inversion parameters calculated by correlating the marginal probability densities of two parameters and normalizing them. Range r and source depth zs correlate with the seabed sound speed b c  , which causes multiple peaks in PPD for cb. When the inversion parameters have multiple peaks in marginal probability densities, many sets computed in parallel are necessary to ensure that the parameters converge to the global optimum. Other parameters do not indicate a strong correlation. When there is strong correlation between inverse parameters, the inverse results will be multiply valued. The seabed sound speed c b can be estimated by seabed speed empirical formula [29]: The inverse parameters were range r, source depth z s , and seabed density ρ b . The lower the number of inverse parameters, the faster the calculating speed. The multi-valuedness can be avoided, and the location result is more accurate.
The estimated values are the MAP values; these results are listed in Table 1.
According to the estimated results, the source range r = 7.67 km and z s = 19.74 m are close to the true values, and the errors are less than 3%. The estimated values of the seabed parameters are also in a good agreement with the true values. The source range r, source depth z s , and seabed density ρ b are well estimated, while the seabed sound speed c b is calculated by Equation (31).

The Validity and Robustness of the Algorithm
To evaluate the effectiveness and robustness of the algorithm, location performance is quantified in terms of the root mean squared error (RMSE) at a variety of SNRs-−10, −5, 0, 5, 10, 15, 20 and 25 dB-in 50 simulation sets. The RMSE can be computed by [31] where M 0 is the number of simulation sets, X i is the inverse results, and X 0 is the true value. In different SNRs, the RMSE of location results (source range and source depth) at a variety of SNRs in 50 simulation sets are shown in Figure 8. Figure 8 shows that the RMSE values decrease with the SNR. When the SNR is higher than 10 dB, the RMSE values near 0. The location results are almost all acceptable. where M0 is the number of simulation sets, Xi is the inverse results, and X0 is the true value. In different SNRs, the RMSE of location results (source range and source depth) at a variety of SNRs in 50 simulation sets are shown in Figure 8. Figure 8 shows that the RMSE values decrease with the SNR. When the SNR is higher than 10 dB, the RMSE values near 0. The location results are almost all acceptable.
(a) The RMES of source range ( b) The RMES of source depth

Processing Results of Measured Data
The experiment was performed in a shallow water waveguide with a half-infinite liquid seabed in an area of the Yellow Sea, China. The depth is 25 m. The sound speed profile in water is shown in Figure 9. A broadband linear frequency modulated impulse source was emitted by UW350 at depth zs = 10 m with a frequency band of 200~500 Hz. The signal duration was 3 s, and the signal was received at depth zr = 9 m by a single sensor after propagation range r = 4972 m, which was measured via a global positioning system (GPS). The seabed sound speed cb and seabed density ρb have been inversed by other methods before; the results will here be compared with these inversed results. The equipment of this experiment is shown in Figure 10.
The pulse compression technique is used to receive signals for convenience. The received signals in the time domain and the TF domain, the warped signal in the TF domain, and the warped signal spectrum are shown in Figures 11a-d, respectively. There are three obvious normal mode signals

Processing Results of Measured Data
The experiment was performed in a shallow water waveguide with a half-infinite liquid seabed in an area of the Yellow Sea, China. The depth is 25 m. The sound speed profile in water is shown in Figure 9. A broadband linear frequency modulated impulse source was emitted by UW350 at depth z s = 10 m with a frequency band of 200~500 Hz. The signal duration was 3 s, and the signal was received at depth z r = 9 m by a single sensor after propagation range r = 4972 m, which was measured via a global positioning system (GPS). The seabed sound speed c b and seabed density ρ b have been inversed by other methods before; the results will here be compared with these inversed results. The equipment of this experiment is shown in Figure 10.   The pulse compression technique is used to receive signals for convenience. The received signals in the time domain and the TF domain, the warped signal in the TF domain, and the warped signal spectrum are shown in Figure 11a-d, respectively. There are three obvious normal mode signals shown in Figure 11b, while four modes can be seen in Figure 11c,d, which illustrate the advantage of the warping transformation. From Figure 11d, we see that the 3rd mode is faint; one reason is that the received sensor is at the null point of the 3rd mode, and the thermocline in the sound speed profile (SSP) may be another reason. To obtain a better estimated result, the 1st, 2nd, and 4th modes with strong energy were used. The extracted modes are shown in Figure 12.
As in Section 4.3, the inversion results are listed in Table 2, and the marginal probability densities of the inversion parameters are shown in Figure 13.  As in Section 4.3, the inversion results are listed in Table 2, and the marginal probability densities of the inversion parameters are shown in Figure 13.   According to the estimated results, the source range  r = 5.04 km was close to the independent range from the GPS data obtained during the experiment (4.792 km), while  s z = 9.99 m was close to the true value, and the errors of the source range and depth were both less than 10%. The estimated values of the seabed parameters were also in a good agreement with the mean inverted values from other methods. The estimated result was acceptable.

Conclusions
This paper presents a passive location method of an underwater source based on a single received sensor. Bayesian methodology was used to build the cost function. The methodology was adapted for low frequency propagation in shallow water. In this case, propagation was dispersive, and the received signal consisted of several normal modes, which could be described in the TF According to the estimated results, the source range r = 5.04 km was close to the independent range from the GPS data obtained during the experiment (4.792 km), while z s = 9.99 m was close to the true value, and the errors of the source range and depth were both less than 10%. The estimated values of the seabed parameters were also in a good agreement with the mean inverted values from other methods. The estimated result was acceptable.

Conclusions
This paper presents a passive location method of an underwater source based on a single received sensor. Bayesian methodology was used to build the cost function. The methodology was adapted for low frequency propagation in shallow water. In this case, propagation was dispersive, and the received signal consisted of several normal modes, which could be described in the TF domain. A signal processing method called the warping transformation was used to separate and extract different numbers of normal modes, which were then used as the input for the Bayesian inversion scheme. A GA was used to search for the optimum solution, and 10 sets were computed in parallel to ensure the parameters converge to the global optimum.
The key point of this paper is that applying different numbers of normal mode signals instead of array (vertical array or horizontal array) data can reduce cost and avoid measurement errors caused by external environmental factors. Additionally, Bayesian methodology is a good inversion method with a rigorous evaluation of data errors and model parameterization, which can realize geoacoustic inversion with estimated parameter uncertainties, so the source range, depth, and other seabed parameters can be estimated without prior knowledge of the seabed information. The estimated results were in a good agreement with true values and estimated values from other inversion methods.
The inversion parameter sensitivity to the cost function and the inter-relationship of parameters was analyzed. The results illustrate the following: The cost function to every single parameter was effective, and source range r and source depth zs correlated with seabed sound speed c b , while other parameters did not indicate a strong correlation. The seabed sound speed c b can be estimated by the seabed speed empirical formula to solve the multi-valuedness caused by strong correlations between inverse parameters. The effectiveness and robustness of the algorithm were quantified in terms of the root mean squared error (RMSE) at a variety of signal-to-noise ratios (SNRs) in 50 simulation sets.
The RMSE values decrease with the SNR. The method described in this paper was applied to the shallow water waveguide with a half-infinite liquid seabed, and the sound speed in water and the water depth were known, so the application to more complex marine environments can be studied in the future.