Particle Filtering for Localization of Broadband Sound Source Using an Ocean-Bottom Seismometer Sensor

Abstract

Passive source localization is a challenging task for one receiver, and the pressure sensor provides relatively simple information. An ocean-bottom seismometer (OBS) sensor placed on the seafloor surface can provide more information—not only pressure information, but also three-axis (x-, y-, and z-axis) velocity information at the seafloor interface. In this paper, an OBS sensor was used to estimate the position of the broadband sound source in a Pekeris shallow water waveguide with elastic bottom. As the dynamics that characterize ocean acoustic applications are inherently nonlinear, non-Gaussian, and non-stationary processes that quickly vary with space and time, sequential Bayesian filtering, such as particle filtering (PF), is able to adapt to these environmental changes. Simulation results show that the PF method with the vertical wave impedance (the ratio of the pressure and vertical particle velocity) in the frequency domain as a measurement vector is not affected by source depth and source spectrum information, making it more tolerant and more robust than that with pressure in positioning. Experimental data results verified the effectiveness of the PF method with the vertical wave impedance for the localization of the explosive source.

1. Introduction

Localization of an underwater sound source is an important practical problem in underwater acoustics. Of all the methods for source localization, matched-field processing (MFP) attracted a huge amount of interest over the years. MFP is a technique combining hydro-acoustic physics and signal processing technology that made important progress in underwater acoustic positioning. However, the matched field method uses a deterministic model, which is plagued by mismatch problems, including environmental mismatch, statistical mismatch, and system mismatch [1]. Sequential filtering provides a suitable framework for estimating and updating unknown parameters of a system as data become available. Moreover, sequential filtering is demonstrated to be a powerful estimation tool, employing prediction from previous estimates and updates stemming from physical and statistical models that relate acoustic measurements to the unknown parameters [2]. An adaptive model-based approach using the state-space formulation was for the first time implemented by Candy et al. [3]. Coupling a nonlinear optimization algorithm with the extended Kalman filter (EKF)-based ocean acoustic model can solve the source localization problem in a complex ocean environment. This approach is capable of solving the mismatch problem to some extent. Furthermore, Candy et al. [4] provided a model-based Bayesian processor to estimate the bearings of moving sources using horizontally towed array data. The processor uses the particle filtering (PF) as the estimator, which has better tracking performance than classical nonlinear filters (e.g., extended/unscented Kalman filters). A particle filtering method for the estimation of such arrivals was presented by Michalopoulou and Jain [5,6]. Furthermore, successful localization with real data was demonstrated using arrival times and corresponding probability density functions (PDFs) extracted via particle filtering [7]. These researches generally used sound pressure signals received by the hydrophone array. In Reference [8], the vertical specific acoustic impedance (the ratio of the complex pressure field and vertical pressure gradient) was used to sense ocean-bottom geo-acoustic properties in the Pekeris environment using a sequential approach. This approach was applied to data collected off the Senegalese coast. Previous studies were generally based on the fluid seabed, without considering the effects of shear waves. However, the seabed in the real marine environment is generally an elastic medium. Therefore, we derive the vertical wave impedance which the ocean-bottom seismometer (OBS) sensor obtained in Pekeris waveguide with an elastic bottom. Then, a passive ranging method for the broadband acoustic source using vertical wave impedance is proposed in this paper based on the particle filtering approach.

The remainder of this paper is organized as follows: in Section 2, the pressure and vertical wave impedance in the frequency domain are constructed in Pekeris waveguide with elastic bottom. In Section 3, the PF framework for estimating the position of the source is introduced. In Section 4, the positioning performances of the pressure and vertical wave impedance as the measurement vector are simulated and analyzed based on the PF method. In Section 5, the presented estimation localization method is employed to process the experimental data. Conclusions are given in Section 6.

2. Theoretical Modeling

In the cylindrical coordinate system, the Pekeris waveguide with elastic bottom, as shown in Figure 1, was built to obtain the pressure and vertical wave impedance received by the OBS sensor. In this section, the pressure field and velocity field expressions are established in normal mode based on wavenumber integration approaches. Thus, in Section 2.1, we briefly introduce the derivation process of potential functions based on wavenumber integration approaches; then, in Section 2.2, the potential function expressions are represented in normal mode form, and then the pressure and the vertical wave impedance at the seafloor interface are obtained.

2.1. Potential Functions

The depth of the uniform water layer is denoted as H, and the density and the sound speed of water are ρ₁ and c₁, respectively. The elastic bottom is assumed to be homogeneous and semi-infinite. The density, the compressional wave, and the shear wave speed of the bottom are constants that are expressed as ρ₂, c_p, and c_s, respectively. As shown in Figure 1, the horizontal axis is range r, the vertical axis is depth z, the plane z = 0 is considered as the sea surface, and the broadband source is positioned at (0, z_s). The spectrum of the source is denoted as S(ω). The OBS sensor is placed on the seafloor interface as a receiver. The potential function of sound field in water is denoted as φ₁, and potential function of sound field in elastic medium is the sum of the scalar potential function φ₂ and vector potential function $\overrightarrow{\psi }$. In the cylindrically symmetric case, only the component along the θ direction exists for $\overrightarrow{\psi }$, which is denoted as ψ referring to Reference [9]. The derivation process of potential functions is detailed in Appendix A.

The expression of ${\phi }_1,{\phi }_2,\psi$ can be written as

$${\phi }_1\left(r,z\right)=\{\begin{array}{c}{\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{2\sin {\beta }_{1}z}{{\beta }_1}\left[\frac{{\beta }_1\cos {\beta }_1\left(H-z_s\right)-ib{\beta }_{2}K\sin {\beta }_1\left(H-z_s\right)}{{\beta }_1\cos {\beta }_{1}H-ib{\beta }_{2}K\sin {\beta }_{1}H}\right]J_0\left(\xi r\right)\xi d\xi ,0\le z<z_s}\\ {\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{2\sin {\beta }_1{z}_s}{{\beta }_1}\left[\frac{{\beta }_1\cos {\beta }_1\left(H-z\right)-ib{\beta }_{2}K\sin {\beta }_1\left(H-z\right)}{{\beta }_1\cos {\beta }_{1}H-ib{\beta }_{2}K\sin {\beta }_{1}H}\right]J_0\left(\xi r\right)\xi d\xi },z_s\le z<H\end{array},$$

(1)

$${\phi }_2=-{\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{2\xi \sigma K}{{\chi }^2}\frac{2b\sin {\beta }_1{z}_s{e}^{i{\beta }_2\left(z-H\right)}}{{\beta }_1\cos {\beta }_{1}H-ib{\beta }_{2}K\sin {\beta }_{1}H}}{\mathrm{J}}_0\left(\xi r\right)\xi d\xi \begin{array}{cc}& z>H\end{array},$$

(2)

$$\psi ={\displaystyle \underset{0}{\overset{\infty }{\int }}\frac{2i\xi {\beta }_{2}K}{{\chi }^2}\frac{2b\sin {\beta }_1{z}_s{e}^{i\gamma \left(z-H\right)}}{{\beta }_1\cos {\beta }_{1}H-ib{\beta }_{2}K\sin {\beta }_{1}H}}{\mathrm{J}}_1\left(\xi r\right)\xi d\xi \begin{array}{cc}& z>H\end{array},$$

(3)

where ${\beta }_1=\sqrt{k_1^2-{\xi }^2}$, ${\beta }_2=\sqrt{k_2^2-{\xi }^2}$, $\gamma =\sqrt{{\chi }^2-{\xi }^2}$, $k_1=\omega /c_1$, $k_2=\omega /c_p$, $\chi =\omega /c_s$, $b={\rho }_1/{\rho }_2$, $\sigma =\left(2{\xi }^2-{\chi }^2\right)/2\xi$,$K={\chi }^4/\left[4{\xi }^2\left({\sigma }^2+{\beta }_2\gamma \right)\right]$, ξ is horizontal wavenumber, and ω is the angular frequency.

2.2. Pressure and Vertical Wave Impedance

According to the formula

$$\begin{array}{c}{\mathrm{J}}_0\left(\xi r\right)=\frac{1}{2}\left[{\mathrm{H}}_0^{\left(1\right)}\left(\xi r\right)+{\mathrm{H}}_0^{\left(2\right)}\left(\xi r\right)\right]\\ {\mathrm{H}}_0^{\left(2\right)}\left(\xi r\right)=-{\mathrm{H}}_0^{\left(1\right)}\left(\xi r{e}^{i\pi }\right)\end{array},$$

the first integral in Equation (1) can be rewritten as

$${\phi }_1\left(r,z\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{\sin {\beta }_{1}z}{{\beta }_1}\left[\frac{{\beta }_1\cos {\beta }_1\left(H-z_s\right)-ib{\beta }_{2}K\sin {\beta }_1\left(H-z_s\right)}{{\beta }_1\cos {\beta }_{1}H-ib{\beta }_{2}K\sin {\beta }_{1}H}\right]{\mathrm{H}}_0^{\left(1\right)}\left(\xi r\right)\xi d\xi },0\le z<z_s,$$

(4)

with ${\mathrm{H}}_0^{\left(1\right)}(\cdot )$ as the first kind of zero-order Hankel function.

By Cauchy’s method of residues, the integral of Equation (4) will be equal to

$$\begin{array}{ll}{\phi }_1\left(r,z\right)& =\mathrm{Res}\left\{{\mathrm{Z}}_1\left(z,\xi \right){\mathrm{H}}_0^{\left(1\right)}\left(\xi r\right)\xi \right\}-{\displaystyle \underset{branch\begin{array}{c}\begin{array}{c}\end{array}line\end{array}}{\int }{\mathrm{Z}}_1\left(z,\xi \right){\mathrm{H}}_0^{\left(1\right)}\left(\xi r\right)\xi d\xi }\\ & ={\phi }_N\left(r,z\right)+{\phi }_L\left(r,z\right)\end{array},$$

(5)

where ${\phi }_N\left(r,z\right)$ is called the normal mode, and ${\phi }_L\left(r,z\right)$ is called the lateral wave.

The expression of ${\phi }_N\left(r,z\right)$ is written as

$$\begin{array}{ll}{\phi }_N\left(r,z\right)& =2\pi i{\displaystyle \sum _n\frac{\sin {\beta }_{1n}z}{{\beta }_{1n}}\frac{{\beta }_{1n}\cos {\beta }_{1n}\left(H-z_s\right)-ib{\beta }_{2n}K\sin {\beta }_{1n}\left(H-z_s\right)}{\frac{\partial }{\partial {\xi }_n}\left({\beta }_{1n}\cos {\beta }_{1n}H-ib{\beta }_{2n}K\sin {\beta }_{1n}H\right)}{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{n}r\right){\xi }_n}\\ & =2\pi i{\displaystyle \sum _n{F}_n^2\sin {\beta }_{1n}z\sin {\beta }_{1n}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{n}r\right)}\end{array}$$

(6)

where

$$F_n^2=\frac{{\beta }_{1n}}{H{\beta }_{1n}-\sin {\beta }_{1n}H\cos {\beta }_{1n}H-b^2{K}_n^2\tan {\beta }_{1n}H{\sin }^2{\beta }_{1n}H-ib{\beta }_{1n}{\beta }_{2n}\frac{\partial K_n}{\partial {\xi }_n}{\sin }^2{\beta }_{1n}H/{\xi }_n}$$

$$\frac{\partial K_n}{\partial {\xi }_n}=-\frac{{\chi }^4}{4}\frac{2{\xi }_n\left({\sigma }_n^2+{\beta }_{2n}{\gamma }_n\right)+{\xi }_n^2\left[2{\sigma }_n\left(1+\frac{{\chi }^2}{2{\xi }_n^2}\right)-{\beta }_{2n}\frac{{\xi }_n}{{\gamma }_n}-\frac{{\xi }_n}{{\beta }_{2n}}{\gamma }_n\right]}{{\xi }_n{}^4{\left({\sigma }_n{}^2+{\beta }_{2n}{\gamma }_n\right)}^2}$$

and values of ${\xi }_n$ are determined by the dispersion equation, i.e., ${\beta }_{1n}\cos {\beta }_{1n}H-ib{\beta }_{2n}{K}_n\sin {\beta }_{1n}H=0$.

Similarly, ${\phi }_N\left(r,z\right)$ in the depth domain $z_s\le z<H$ also can be obtained with the same expression as Equation (6).

Due to the effect of the shear wave, the expression of the lateral wave is more complicated than that in Reference [9]. Specifically, the lateral wave mainly consists of two parts. One integration along the branch at first runs from $k_2+i\infty$ to $k_2$ along the left side of the branch line with the negative imaginary part of ${\beta }_2$, then runs from $k_2$ to $k_2+i\infty$ along the right side of the branch line with the positive imaginary part of ${\beta }_2$. The other integration along the branch at first runs from $\chi +i\infty$ to $\chi$ along the left side of the branch line with the negative imaginary part of $\gamma$, then runs from $\chi$ to $\chi +i\infty$ along the right side of the branch line with the positive imaginary part of $\gamma$. Under the assumption that $k_2<\chi <k_1$, the approximate expression of lateral wave is

$${\phi }_L\left(r,z\right)={\phi }_{Lk}\left(r,z\right)+{\phi }_{L\chi }\left(r,z\right),$$

(7)

where

$$\begin{array}{ll}{\phi }_{Lk}\left(r,z\right)& ={\displaystyle {\int }_{k_2}^{k_2+i\infty }\frac{\sin {\beta }_{1}z}{{\beta }_1}}\{\frac{{\left(2{\xi }^2-{\chi }^2\right)}^2{\beta }_1\cos {\beta }_1\left(H-z_s\right)+{\beta }_2\left[4{\xi }^2\gamma {\beta }_1\cos {\beta }_1\left(H-z_s\right)-ib{\chi }^4\sin {\beta }_1\left(H-z_s\right)\right]}{{\left(2{\xi }^2-{\chi }^2\right)}^2{\beta }_1\cos {\beta }_{1}H+{\beta }_2\left[4{\xi }^2\gamma {\beta }_1\cos {\beta }_{1}H-ib{\chi }^4\sin {\beta }_{1}H\right]}\\ & -\frac{{\left(2{\xi }^2-{\chi }^2\right)}^2{\beta }_1\cos {\beta }_1\left(H-z_s\right)-{\beta }_2\left[4{\xi }^2\gamma {\beta }_1\cos {\beta }_1\left(H-z_s\right)-ib{\chi }^4\sin {\beta }_1\left(H-z_s\right)\right]}{{\left(2{\xi }^2-{\chi }^2\right)}^2{\beta }_1\cos {\beta }_{1}H-{\beta }_2\left[4{\xi }^2\gamma {\beta }_1\cos {\beta }_{1}H-ib{\chi }^4\sin {\beta }_{1}H\right]}\}{\mathrm{H}}_0^{\left(1\right)}\left(\xi r\right)\xi d\xi \end{array}$$

$$\approx \{\begin{array}{c}\begin{array}{ccc}& & \end{array}\frac{2b{k}_2{\chi }^4\sin \sqrt{k_1^2-k_2^2}{z}_0\sin \sqrt{k_1^2-k_2^2}z}{\left(k_1^2-k_2^2\right){\cos }^2\sqrt{k_1^2-k_2^2}H}\frac{1}{r^2}{e}^{i\left(k_{2}r-\frac{\pi }{2}\right)},\begin{array}{cc}& \begin{array}{cc}& \end{array}\end{array}if\cos \sqrt{k_1^2-k_2^2}H\ne 0\\ -\frac{2b{\left(2k_2{}^2-{\chi }^2\right)}^2{\chi }^4\sin \sqrt{k_1^2-k_2^2}{z}_0\sin \sqrt{k_1^2-k_2^2}z}{\left[4k_2^2\gamma \sqrt{k_1^2-k_2^2}\cos \sqrt{k_1^2-k_2^2}H-ib{\chi }^4\sin \sqrt{k_1^2-k_2^2}H\right]}\frac{1}{r}{e}^{i{k}_{2}r},\begin{array}{cc}& \end{array}if\cos \sqrt{k_1^2-k_2^2}H=0\end{array},$$

(8)

$$\begin{array}{l}{\phi }_{L\chi }\left(r,z\right)={\displaystyle {\int }_{\chi }^{\chi +i\infty }\frac{\sin {\beta }_{1}z}{{\beta }_1}\{\frac{{\left(2{\xi }^2-{\chi }^2\right)}^2{\beta }_1\cos {\beta }_1\left(H-z_s\right)-ib{\beta }_2{\chi }^4\sin {\beta }_1\left(H-z_s\right)+4{\xi }^2{\beta }_2\gamma {\beta }_1\cos {\beta }_1\left(H-z_s\right)}{{\left(2{\xi }^2-{\chi }^2\right)}^2{\beta }_1\cos {\beta }_{1}H-ib{\beta }_2{\chi }^4\sin {\beta }_{1}H+4{\xi }^2{\beta }_2\gamma {\beta }_1\cos {\beta }_{1}H}}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{{\left(2{\xi }^2-{\chi }^2\right)}^2{\beta }_1\cos {\beta }_1\left(H-z_s\right)-ib{\beta }_2{\chi }^4\sin {\beta }_1\left(H-z_s\right)-4{\xi }^2{\beta }_2\gamma {\beta }_1\cos {\beta }_1\left(H-z_s\right)}{{\left(2{\xi }^2-{\chi }^2\right)}^2{\beta }_1\cos {\beta }_{1}H-ib{\beta }_2{\chi }^4\sin {\beta }_{1}H-4{\xi }^2{\beta }_2\gamma {\beta }_1\cos {\beta }_{1}H}\}{\mathrm{H}}_0^{\left(1\right)}\left(\xi r\right)\xi d\xi \\ \approx \{\begin{array}{c}\frac{8\sin \sqrt{k_1^2-{\chi }^2}{z}_s\sin \sqrt{k_1^2-{\chi }^2}z}{b\chi {\sin }^2\sqrt{k_1^2-{\chi }^2}H}\frac{1}{r^2}{e}^{i\left(\chi r+\frac{\pi }{2}\right)},if\cos \sqrt{k_1^2-{\chi }^2}H=0\\ -\frac{{\chi }^{2}b\sin \sqrt{k_1^2-{\chi }^2}{z}_s\sin \sqrt{k_1^2-{\chi }^2}z}{2\left(k_1^2-{\chi }^2\right){\cos }^2\sqrt{k_1^2-{\chi }^2}H}\frac{1}{r}{e}^{i\chi r},if\left(\begin{array}{c}\cos \sqrt{k_1^2-{\chi }^2}H\ne 0,\\ b\sqrt{{\chi }^2-k_2^2}\tan \sqrt{k_1^2-{\chi }^2}H+\sqrt{k_1^2-{\chi }^2}=0\end{array}\right)\\ \frac{8\left(k_2^2-{\chi }^2\right)b\sin \sqrt{k_1^2-{\chi }^2}{z}_s\sin \sqrt{k_1^2-{\chi }^2}z}{\chi {\left[\sqrt{k_1^2-{\chi }^2}\cos \sqrt{k_1^2-{\chi }^2}H+b\sqrt{{\chi }^2-k_2^2}\sin \sqrt{k_1^2-{\chi }^2}H\right]}^2}\frac{1}{r^2}{e}^{i\left(\chi r-\frac{\pi }{2}\right)},if\left(\begin{array}{c}\cos \sqrt{k_1^2-{\chi }^2}H\ne 0,\\ b\sqrt{{\chi }^2-k_2^2}\tan \sqrt{k_1^2-{\chi }^2}H+\sqrt{k_1^2-{\chi }^2}\ne 0\end{array}\right)\end{array}\end{array}$$

(9)

Equations (8) and (9) show that, in general, lateral waves attenuate more quickly than spherical waves during propagation, and will possess a significant value only for distances not far from the source. In order to reflect the phenomenon more intuitively, the transmission losses (TLs) of normal modes, the lateral wave, and the spherical wave are shown in Figure 2. The simulation environment parameters are shown in

Table 1. Parameters of the ocean environment.

Medium	Depth(m)	Density(g/cm3)	Compression Wave Speed (m/s)	Shear Wave Speed (m/s)
Fluid	50	1.0	1500	/
Elastic bottom	/	1.5	3800	1800

, the sound source depth is 20 m, and the receiver is placed on the seafloor interface. Here, the frequency is 50 Hz. Figure 2 shows that the lateral wave attenuation decays faster with distance. Thus, in general, when long-range sound transmission is considered, the contribution of the lateral wave usually can be neglected and

$${\phi }_1\left(r,z\right)\approx {\phi }_N\left(r,z\right)=2\pi i{\displaystyle \sum _n{F}_n^2\sin {\beta }_{1n}z\sin {\beta }_{1n}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{n}r\right)}.$$

(10)

In the fluid, the relationships between potential function, complex pressure, and the vertical particle velocity at angular frequency ω are as follows:

$$p\left(r,z,\omega \right)={\rho }_1{\omega }^{2}S\left(\omega \right){\phi }_1\left(r,z\right)=2\pi iS\left(\omega \right){\rho }_1{\omega }^2{\displaystyle \sum _n{F}_n^2\sin {\beta }_{1n}z\sin {\beta }_{1n}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{n}r\right)};$$

(11)

$$v_z\left(r,z,\omega \right)=-i\omega S\left(\omega \right)\frac{\partial {\phi }_1\left(r,z\right)}{\partial z}=\omega 2\pi S\left(\omega \right){\displaystyle \sum _n{F}_n^2{\beta }_{1n}\cos {\beta }_{1n}z\sin {\beta }_{1n}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{n}r\right)}$$

(12)

Because the vertical particle velocity is continuous at z = H, the ratio of the pressure and the vertical velocity received by the OBS is

$$Z_z\left(r,H,\omega \right)=\frac{p\left(r,H,\omega \right)}{v_z\left(r,H,\omega \right)}=i{\rho }_1\omega \frac{{\displaystyle \sum _n{F}_n^2\sin {\beta }_{1n}H\sin {\beta }_{1n}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{n}r\right)}}{{\displaystyle \sum _n{F}_n^2{\beta }_{1n}\cos {\beta }_{1n}H\sin {\beta }_{1n}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{n}r\right)}}.$$

(13)

If only the single normal mode is considered, the expression of vertical wave impedance is represented as (Z_z = iρ₁ωsinβ_1nH/β_1ncosβ_1nH). Hence, the vertical wave impedance is only related to the receiver depth and is independent of the distance. Therefore, for ranging, the vertical wave impedance must take into account the normal mode order of two and above.

3. PF Framework

PF is a sequential Monte Carlo (MC) method employing the sequential estimation of relevant probability distributions using the importance sampling (IS) techniques and the approximations of distributions with discrete random measures [10,11,12,13]. PF is a technique to implement sequential Bayesian estimators via MC simulation [14]. PF has the advantage that the noise model and the system model are not limited. Compared with the Kalman filter, PF can be applied to nonlinear systems and is not limited by Gaussian distribution. The PF framework for estimating the position of the source is organized as follows: in Section 3.1, a general background about the state-space equation for the estimation of evolving parameters in dynamical systems is given together with the basics of Bayesian filtering. In Section 3.2, the filter equations are derived starting from the basic IS concepts, moving to sequential importance sampling (SIS), finally deriving the commonly used PF, often referred to as sequential importance resampling (SIR).

3.1. Background

3.1.1. State-Space Model

The state-space model requires two equations: the state equation and the measurement equation. The state equation models the evolution of the state vector over time. The measurement equation performs the mapping from state vectors to observation vectors. The two equations are represented as follows:

$${\mathbf{x}}_k=\mathbf{f}\left({\mathbf{x}}_{k-1}\right)+{\mathbf{w}}_{k-1},$$

(14)

$${\mathbf{y}}_k=\mathbf{d}\left({\mathbf{x}}_k\right)+{\mathbf{v}}_k,$$

(15)

where x = [r, z_s]^T, r is the range, z_s is the source depth, and T represents the transpose of matrix x.

The state Equation (14) describes the evolution or transition of x_k and assumes that states follow a first-order Markov process. Function f(∙) is the state prediction operator, which models the evolution of the state vector at step k − 1 to that of step k. In this work, f(∙) is taken as the unit matrix I. The measurement Equation (15) relates measurements y_k to state vector x_k through function d(∙), where y is the measurement vector. Function d(∙) is the nonlinear function that relates the environmental and source parameters x_k to the acoustic measurement vector y_k. In this work, when vertical wave impedance selected as measurement vector, y = [Z_z(x, ω₁), …, Z_z (x, ω_m), …, Z_z (x, ω_M)]^T, and Z_z (x, ω_m) is the vertical wave impedance at angular frequency ω_m, with $m\in \left[1,\mathrm{M}\right]$. Furthermore, d(x) = [d₁(x), …, d_m(x), …, d_M(x)]^T, with

$${\mathbf{d}}_m\left(\mathbf{x}\right)=i{\rho }_1{\omega }_m\frac{{\displaystyle \sum _n{F}_{nm}^2\sin {\beta }_{1nm}H\sin {\beta }_{1nm}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{nm}r\right)}}{{\displaystyle \sum _n{F}_{nm}^2{\beta }_{1nm}\cos {\beta }_{1nm}H\sin {\beta }_{1nm}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{nm}r\right)}}.$$

(16)

Moreover, w_k and v_k are the state noise vector and the measurement noise vector, respectively, with

$$E\left\{{\mathbf{w}}_k{\mathbf{w}}_i^T\right\}={\mathbf{Q}}_k{\delta }_{ki},E\left\{{\mathbf{v}}_k{\mathbf{v}}_i^T\right\}={\mathbf{R}}_k{\delta }_{ki},E\left\{{\mathbf{w}}_k{\mathbf{v}}_i^T\right\}=0,\forall i,k,$$

(17)

where δ(∙) is the Dirac delta function, and Q_k and R_k are the covariance matrices at k for the corresponding noise terms. The modal/range uncertainty can be lumped as a state process noise term to represent sound-speed profile errors, errors in the boundary conditions, sea state effects, and ocean inhomogeneities, whereas the measurement noise can be lumped into an additive noise term to represent the near-field acoustic noise field, flow noise on the sensors, and electronic noise [15].

3.1.2. Bayesian Filtering

After giving the state equation and the measurement equation, we briefly discuss the algorithm for our estimation problem. The basic problem we pursue in this paper can now be defined in terms of our mathematical models as follows:

GIVEN a set of noisy vertical wave impedances along with the state equation and measurement equation (Equations (14) and (15)) with unknown parameters (r, z_s), FIND the “best” estimated values.

Examining the problem from a Bayesian standpoint [16,17], we are interested in deriving the full posterior PDF for x_k. The initial PDF of the state vector, Pr(x₀), is assumed to be known. Let D_k = [y₁, y₂, …, y_k] be the set of data from 1 to k steps. The aim is to estimate Pr(x_k|D_k), the posterior PDF of the state vector at step k.

With the posterior PDF Pr(x_{k − 1}|D_{k − 1}) available, we can predict Pr(x_k|D_{k − 1}) through the transition PDF Pr(x_k|x_{k − 1}). Due to the first-order Markov chain assumption of x_k, Pr(x_k|x_{k − 1}) does not depend on data D_{k − 1}. Density Pr(x_k|D_{k − 1}) can be written as

$$\begin{array}{ll}\mathrm{Pr}\left({\mathbf{x}}_k|{\mathbf{D}}_{k-1}\right)& ={\displaystyle \int \mathrm{Pr}\left({\mathbf{x}}_k|{\mathbf{x}}_{k-1},{\mathbf{D}}_{k-1}\right)}\times \mathrm{Pr}\left({\mathbf{x}}_{k-1}|{\mathbf{D}}_{k-1}\right)d{\mathbf{x}}_{k-1}\\ & ={\displaystyle \int \mathrm{Pr}\left({\mathbf{x}}_k|{\mathbf{x}}_{k-1}\right)}\times \mathrm{Pr}\left({\mathbf{x}}_{k-1}|{\mathbf{D}}_{k-1}\right)d{\mathbf{x}}_{k-1}\end{array}$$

(18)

When a new measurement y_k becomes available, the posterior PDF Pr(x_k|D_k) can be calculated by the Bayes theorem.

$$\begin{array}{ll}\mathrm{Pr}\left({\mathbf{x}}_k|{\mathbf{D}}_k\right)& =\frac{\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)\mathrm{Pr}\left({\mathbf{x}}_{k-1}|{\mathbf{D}}_{k-1}\right)}{\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{D}}_{k-1}\right)}\\ & =\frac{\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)\mathrm{Pr}\left({\mathbf{x}}_{k-1}|{\mathbf{D}}_{k-1}\right)}{{\displaystyle \int \mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)\mathrm{Pr}\left({\mathbf{x}}_{k-1}|{\mathbf{D}}_{k-1}\right)}d{\mathbf{x}}_{\mathrm{k}}}\end{array}$$

(19)

The posterior PDF Pr(x_k|D_k) contains all information provided from the data, the measurement equation, and the noise model about target x_k at step k.

3.2. Particle Filter

PFs track the posterior PDF Pr(x_k|D_k) using a cloud of particles ${\left\{{\mathbf{x}}_k^i\right\}}_{i=1}^{N_p}=\left\{{\mathbf{x}}_k^1,{\mathbf{x}}_k^2,\cdots ,{\mathbf{x}}_k^{N_p}\right\}$ that evolve with step k. Before presenting the details of the PF, we summarize the basic IS concepts.

3.2.1. Importance Sampling

For a general nonlinear system, it is difficult to obtain an analytical solution of the posterior probability, and it is difficult to obtain the integral in Equations (18) and (19). In order to solve the integral problem, we introduce the MC method.

Assume that we want to compute an integral I = ∫f(x)dx. One way of computing I is assuming x is a random variable, defining f(x) = g(x)p(x), and rewriting it in the form of an expectation [2].

$$I=E\left\{g\left(x\right)\right\}={\displaystyle \int g}\left(x\right)p\left(x\right)dx,$$

(20)

where g(x) is some function of x with PDF p(x). By drawing N_p independent and identically distributed x samples from p(x), I can be computed numerically via MC integration [2].

$${\left\{x^i\right\}}_{i=1}^{N_p}\sim p\left(x\right)\to I\approx \frac{1}{N_p}{\displaystyle \sum _{i=1}^{N_p}g\left(x^i\right)}.$$

(21)

However, in many cases, it is too costly or not possible to sample from p(x). To mitigate difficulties with inability to directly from a posterior distribution, IS is introduced. IS is a method to compute expectations with respect to one density using random samples drawn from another. Using a simple function q(x) as the sampling density, Equation (20) can be rewritten as

$$\begin{array}{l}I={\displaystyle \int \left[g\left(x\right)\frac{p\left(x\right)}{q\left(x\right)}\right]}q\left(x\right)dx=E\left\{g\left(x\right)\frac{p\left(x\right)}{q\left(x\right)}\right\}\\ \begin{array}{cccc}& & & \end{array}{\left\{x^i\right\}}_{i=1}^{N_p}\sim q\left(x\right)\end{array}$$

(22)

The estimate is obtained using MC integration,

$$\widehat{I}=\frac{1}{N_p}{\displaystyle \sum _{i=1}^{N_p}g\left(x^i\right)\frac{p\left(x^i\right)}{q\left(x^i\right)}}=\frac{1}{N_p}{\displaystyle \sum _{i=1}^{N_p}{w}^{i}g\left(x^i\right)},$$

(23)

where $w^i=p\left(x^i\right)/q\left(x^i\right)$ represents the importance weights.

3.2.2. Sequential Importance Sampling

Bayesian filtering requires performing successive IS runs at each k. The output of each IS run is used as the prior for the next one. This process is referred to as SIS. Let X_k = [x₁, x₂, …, x_k]; it is possible to obtain posterior PDF Pr(x_k|D_k) from the full posterior density Pr(X_k|D_k).

$$\mathrm{Pr}\left({\mathbf{x}}_k|{\mathbf{D}}_k\right)={\displaystyle \int \delta \left({\mathbf{x}}_k-{\mathbf{x}}_k{}^{\prime }\right)}\mathrm{Pr}\left(X_k{}^{\prime }|{\mathbf{D}}_k\right)\mathrm{d}{X}_k{}^{\prime }.$$

(24)

Selecting a sampling density q(X_k|D_k) and implementing IS, we can obtain

$$\mathrm{Pr}\left({\mathbf{x}}_k|{\mathbf{D}}_k\right)\approx {\displaystyle \sum _{i=1}^{N_p}{W}_k^i\delta \left({\mathbf{x}}_k-{\mathbf{x}}_k^i\right)},$$

(25)

$$W_k^i\propto \frac{\mathrm{Pr}\left({\mathbf{X}}_k^i|{\mathbf{D}}_k^i\right)}{q\left({\mathbf{X}}_k^i|{\mathbf{D}}_k^i\right)}.$$

(26)

Expanding the full posterior PDF [10], Pr(X_k|D_k) can be expressed as

$$\mathrm{Pr}\left({\mathbf{X}}_k|{\mathbf{D}}_k\right)=\frac{\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k\right)\mathrm{Pr}\left({\mathbf{x}}_k|{\mathbf{x}}_{k-1}\right)}{\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{D}}_{k-1}\right)}\mathrm{Pr}\left({\mathbf{X}}_{k-1}|{\mathbf{D}}_{k-1}\right).$$

(27)

Selecting q(X_k|D_k) as

$$q\left({\mathbf{X}}_k|{\mathbf{D}}_k\right)=q\left({\mathbf{x}}_k|{\mathbf{x}}_{k-1},{\mathbf{D}}_k\right)q\left({\mathbf{X}}_{k-1}|{\mathbf{D}}_{k-1}\right),$$

(28)

the weight of the ith particle at step k can be represented as

$$\begin{array}{ll}{W}_k^i& \propto \frac{\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k^i\right)\mathrm{Pr}\left({\mathbf{x}}_k^i|{\mathbf{x}}_{k-1}^i\right)}{q\left({\mathbf{x}}_k^i|{\mathbf{x}}_{k-1}^i,{\mathbf{D}}_k\right)}\times \frac{\mathrm{Pr}\left({\mathrm{X}}_{k-1}|{\mathbf{D}}_{k-1}\right)}{q\left({\mathrm{X}}_{k-1}^i|{\mathbf{D}}_{k-1}^i\right)}\\ & \propto \frac{\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k^i\right)\mathrm{Pr}\left({\mathbf{x}}_k^i|{\mathbf{x}}_{k-1}^i\right)}{q\left({\mathbf{x}}_k^i|{\mathbf{x}}_{k-1}^i,{\mathbf{D}}_k\right)}{W}_{k-1}^i\end{array}$$

(29)

There are a variety of the PF algorithms available, each evolving from a particular choice of sampling density; however, perhaps the simplest is the bootstrap technique [18], which we apply to our problem. Here, the sampling density is selected as

$$q\left({\mathbf{x}}_k^{}|{\mathbf{x}}_{k-1}^{},{\mathbf{D}}_k\right)=\mathrm{Pr}\left({\mathbf{x}}_k^{}|{\mathbf{x}}_{k-1}^{}\right).$$

(30)

Thus, Equation (29) is simplified to

$$W_k^i=\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k^i\right)W_{k-1}^i.$$

(31)

Now, only Pr(y_k|x_k) at step k is employed in updating weights $W_k^i$. Pr(y_k|x_k) is called the likelihood function.

3.2.3. Sequential Importance Resampling

One of the major problems with SIS is the degeneracy of the particles. After a few iterations of successive SIS, the process leads to a cloud containing few particles with large weights and numerous particles with negligible weights. This loss of sample diversity results in poor filter performance. Thus, there is a need to somehow resolve this degeneracy problem. This requirement leads to the idea of “resampling” the particles. SIS with an additional resampling step to avoid degeneracy called SIR [16,17,19], and SIR is the most popular PF implementation.

Resampling is easily performed at the end of each step k; alternatively, resampling is implemented when the effective number of particles N_eff needed to maintain diversity drops below a threshold N_thresh. An estimate of the effective number of particles is given by

$$N_{eff}=\frac{1}{{\displaystyle \sum _{i=1}^{N_p}{\left(W_k^i\right)}^2}}.$$

(32)

When N_eff is less than the threshold, resampling is performed.

$$N_{eff}=\{\begin{array}{c}\le N_{thresh}\begin{array}{c}\end{array}\mathrm{Re}\mathrm{sample}\\ >N_{thresh}\begin{array}{c}\end{array}\begin{array}{c}\end{array}\mathrm{Accept}\end{array}.$$

(33)

In summary, the SIR particle filter works as follows: suppose that, at time k – 1, there is a particle cloud $\left\{{\mathbf{x}}_{k-1}^1,{\mathbf{x}}_{k-1}^2,\cdots ,{\mathbf{x}}_{k-1}^{N_p}\right\}$ of size N_p that, with associated weights, samples from the posterior PDF Pr(x_{k − 1}|D_{k − 1}). Then, transform cloud $\left\{{\mathbf{x}}_{k-1}^1,{\mathbf{x}}_{k-1}^2,\cdots ,{\mathbf{x}}_{k-1}^{N_p}\right\}$ into $\left\{{\mathbf{x}}_k^{1\ast },{\mathbf{x}}_k^{2\ast },\cdots ,{\mathbf{x}}_k^{N_p\ast }\right\}$ through the state equation. Each article in the latter cloud has weight 1/N_p. When data y_k are available, the normalized weight of each particle is reevaluated, i.e., $W_k^i=\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k^{i\ast }\right)/{\displaystyle {\sum }_{j=1}^{N_p}\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k^{j\ast }\right)}$, where PDF $\mathrm{Pr}\left({\mathbf{y}}_k|{\mathbf{x}}_k^{j\ast }\right)$ is defined by the measurement equation and knowledge of the statistical behavior of errors in data measurements. These weights are used for the estimation of posterior PDF Pr(x_k|D_k) [20].

Thus, once the posterior is available, the estimates of important statistics can be inferred. For instance, the minimum mean-squared error (MMSM) estimate is used in this paper, with

$${\hat{\mathrm{x}}}_{MMSE}\approx \frac{1}{N_p}{\displaystyle \sum _{i=1}^{N_p}{W}_i{\mathbf{x}}_i}.$$

(34)

4. Simulation

The simulation was performed in a shallow water waveguide with a half-infinite elastic seabed as shown in Figure 1. The simulation environment parameters are shown in Table 1. The receiver was placed on the seafloor interface, and the range was 10 km. For convenience, S(ω) was constant.

Firstly, the effects of the particle size and the resampling strategies on the estimation performance were analyzed by simulation. Here, three widely used resampling algorithms (multinomial resampling, systematic resampling, and residual resampling) are discussed when the numbers of PF particles were 150, 300, 600, and 1200. The sound source depth was 20 m, the vertical wave impedance was the measured vector, and the frequency domain was selected as 50 Hz–150 Hz. The range estimation results and source depth estimation results in different resampling algorithms and particle sizes are shown in

Table 2. Range estimation results in different resampling algorithms and particle sizes.

Particle Size	Multinomial Resampling	Systematic Resampling	Residual Resampling
150	9904.4 m	9869.6 m	9902.5 m
300	9904.0 m	9910.4 m	9894.5 m
600	9957.9 m	9907.1 m	9941.0 m
1200	9938.0 m	9949.1 m	9950.3 m

and

Table 3. Source depth estimation results in different resampling algorithms and particle sizes.

Particle Size	Multinomial Resampling	Systematic Resampling	Residual Resampling
150	25.6 m	28.8 m	25.4 m
300	25.7 m	25.2 m	26.3 m
600	22.6 m	26.3 m	23.7 m
1200	22.6 m	22.1 m	21.1 m

, respectively. It should be noted that each of the estimation values in Table 2 and Table 3 is the average of 20 simulation results. It can be seen from Table 2 and Table 3 that, as the number of particles increases, the positioning performances of the three resampling algorithms tend to increase. Following comprehensive comparison of the positioning results of the three resampling algorithms, the best one was determined as residual resampling. Thus, in all subsequent simulations and experiments, if not specified, the resampling algorithm selected residual resampling and the number of particles was chosen to be 1200.

Next, in the same simulation environment, the effect of Q_k (the covariance of state process noise) on the estimation performance was analyzed by simulation. As the initial value of x will affect the value of Q_k, in two different initial values, the influence of Q_k on the positioning performance was analyzed and the empirical rule for determining Q_k was obtained through simulation. The two different initial values of x were [9500, 25]^T and [7500, 25]^T.

When the initial value was [9500, 25]^T, the different values of ${\mathit{Q}}_k^{1/2}$ are given in

Table 4. Different values of ${\mathit{Q}}_k^{1/2}$ in different initial values of x.

The Initial Value ofxWas [9500, 25]T
No.	(1a)	(1b)	(1c)	(1d)	(1e)
${\mathit{Q}}_k^{1/2}$	$\left[\begin{array}{cc}{10}^2& 0\\ 0& 0.5\end{array}\right]$	$\left[\begin{array}{cc}{10}^2& 0\\ 0& 1\end{array}\right]$	$\left[\begin{array}{cc}{10}^2& 0\\ 0& 4\end{array}\right]$	$\left[\begin{array}{cc}10& 0\\ 0& 1\end{array}\right]$	$\left[\begin{array}{cc}{10}^3& 0\\ 0& 1\end{array}\right]$
The Initial Value of x Was [7500, 25]T
No.	(2a)	(2b)	(2c)	(2d)	(2e)
${\mathit{Q}}_k^{1/2}$	$\left[\begin{array}{cc}{10}^3& 0\\ 0& 0.5\end{array}\right]$	$\left[\begin{array}{cc}{10}^3& 0\\ 0& 1\end{array}\right]$	$\left[\begin{array}{cc}{10}^3& 0\\ 0& 4\end{array}\right]$	$\left[\begin{array}{cc}{10}^2& 0\\ 0& 1\end{array}\right]$	$\left[\begin{array}{cc}2\ast {10}^3& 0\\ 0& 1\end{array}\right]$

and the localization results are shown in Figure 3. Combined with Figure 3 and Table 4, the following can be determined:

(i) Figure 3a,b show the positioning results when the second values on the diagonal of ${\mathit{Q}}_k^{1/2}$ were 0.5, 1, and 4, and the first value on the diagonal remained 100, which corresponds to cases (1a), (1b), and (1c), respectively. The distance and source depth cannot be correctly estimated in case (1a), but in cases (1b) and (1c), source position can be correctly estimated. At the same time, the source depth estimation curve in case 1(b) had a relatively smaller fluctuation around the true value compared to case 1(c). In all cases of (1a), (1b), and (1c), the positioning performance was best in case (1b);

(ii) Figure 3c,d show the positioning results when the first values on the diagonal of ${\mathit{Q}}_k^{1/2}$ were 10, 10², and 10³, and the second value on the diagonal remained 1, which corresponds to cases (1d), (1b), and (1e), respectively. The distance and source depth cannot be correctly estimated in case (1d), but in cases (1b) and (1e), source position can be correctly estimated. Furthermore, the source depth estimation curve in case 1(e) had a relatively larger fluctuation around the true value compared to case 1(b). In all cases of (1d), (1b), and (1e), the positioning performance was best in case (1b);

(iii) The value of Q_k affects the positioning performance. Covariance matrix Q_k should contain values large enough to accommodate the unexpected changes, but at the same time, the value of covariance matrix should not be too small which can cause poor positioning performance (such as in cases (1a) and (1d)). When the absolute difference between the initial value and true value of x was [500, 5]^T, the positioning performance was best when ${\mathit{Q}}_k^{1/2}$ was $\left[\begin{array}{cc}{10}^2& 0\\ 0& 1\end{array}\right]$ among the five cases.

When the initial value was [7500, 25]^T, the different values of ${\mathit{Q}}_k^{1/2}$ are also given in Table 4 and the localization results are shown in Figure 4. Although values of ${\mathit{Q}}_k^{1/2}$ are different from that in the case of the initial value being [9500, 25]^T, the empirical rule of selecting Q_k obtained in Figure 4 is similar to that in Figure 3. The value of Q_k relates to the absolute difference between the initial value and the true value of x, and the value of Q_k cannot be selected too small, which will lead to failure in source localization as shown in the curves of (2a) and (2d) in Figure 4. Although the distance and source depth can be correctly estimated when the value of Q_k is large enough, the estimated curves may have relatively larger fluctuations around the true value, for example, in the case of ${\mathit{Q}}_k^{1/2}$ being $\left[\begin{array}{cc}2\ast {10}^3& 0\\ 0& 1\end{array}\right]$. In all cases of (2a)–(2e), the positioning performance was best when ${\mathit{Q}}_k^{1/2}$ was $\left[\begin{array}{cc}{10}^3& 0\\ 0& 1\end{array}\right]$, i.e., case (2b). There was no fixed value for the covariance matrix Q_k, but in the following simulations and experiments, the Q_k value could be determined according to the empirical rule obtained from these simulations.

Then, in the same simulation conditions, the estimation performances of EKF and the unscented Kalman filter (UKF) were compared with PF. The localization results are shown in Figure 5. It can be seen from Figure 5a that UKF and PF performed better than EKF in terms of range estimation. However, in terms of source depth estimation, both UKF and EKF estimation curves had larger fluctuations near the true value (20 m) than the PF estimation curve, as shown in Figure 5b. Figure 5 shows that the positioning performance of PF is superior to both EKF and UKF.

Then, under different source depths conditions, we discussed the source localization performance of pressure and vertical wave impedance as the measurement vector.

4.1. Source Depth of 20 m

When the sound source depth was 20 m, the measured vectors were the pressure and vertical wave impedance. The frequency domain was selected as 25 Hz–100 Hz, and the localization results are shown in Figure 6. For both pressure and vertical wave impedance, the range estimation results converged and were in good agreement with the true range of 10 km (Figure 6a,c). As can be seen from Figure 6b and Figure 6d, both depth iteration curves converged to the true value of 20 m; however, compared to Figure 6d, the curve in Figure 6b (the pressure as the measurement vector) had large fluctuation near the true source depth. In

Table 5. Mean absolute percentage error (MAPE) of localization results in the frequency domain of 25 Hz–100 Hz.

	Range	Source Depth
Pressure	0.0023%	2.7709%
Vertical wave impedance	0.0107%	0.5129%

, the mean absolute percentage error (MAPE) values of localization results are given for comparing the positioning performance between pressure and vertical wave impedance. MAPE was calculated as the average of the unsigned percentage error, as shown in the example below.

$$\mathrm{MAPE}=\left(\frac{1}{N}{\displaystyle \sum _{n=1}^N\frac{\left|{\widehat{x}}_n-x_{true}\right|}{\left|x_{true}\right|}}\right)\ast 100,$$

(35)

where N is the number of estimates. Here, the estimate values were selected from iterations 20 to 100. It can be clearly seen from Table 5 that the pressure was slightly better than the vertical wave impedance when estimating the distance, but neither MAPE exceeded 0.02%. Conversely, the vertical wave impedance was better in estimating the source depth than pressure.

Next, different frequency bands (50 Hz–150 Hz) were selected for the simulation, and the localization results are shown in Figure 7. The MAPE values of localization results are given in

Table 6. MAPE of localization results in the frequency domain of 50 Hz–150 Hz.

	Range	Source Depth
Pressure	0.0031%	2.2607%
Vertical wave impedance	0.0117%	0.6860%

. At 50 Hz–150 Hz, the positioning performances of pressure and vertical wave impedance were basically the same as those at 25 Hz–100 Hz except that the MAPE values were slightly different. Under the simulation conditions, Figure 6 and Figure 7, and Table 5 and Table 6 prove that the positioning performances of pressure and vertical wave impedance were not affected by the frequency band too much.

4.2. Source Depth of 40 m

When the sound source depth was 40m, the measured vectors were the pressure and vertical wave impedance. The frequency domain was selected as 75 Hz–95 Hz, and the localization results are shown in Figure 8. In this case, the performance of vertical wave impedance for localization was still good, as shown in Figure 8c,d. However, when the pressure was used for distance estimation, the result was extremely poor, as shown in Figure 8a, where the iteration estimation curve had great fluctuation at the true value of 10 km. Figure 8b shows that the source depth estimation results converged to the true value of 40 m, but the curve fluctuated within 35 m–45 m. Upon increasing the number of iterations and simulating in the same environment, the localization results of the pressure are shown in Figure 9. Increasing the number of iterations still could not change the ranging capability of the pressure.

The reason for the appearance of Figure 8 is given from the perspective of sensitivity analysis. The effect of source localization depends on the sensitivity of parameters. The sensitivity of parameters refers to the change degree of the objective function caused by the change of the parameter. The normalized objective function used here is given by

$$\mathsf{\Phi }\left(\mathbf{x}\right)=1-\frac{1}{{\mathbf{y}}^{\mathrm{H}}\mathbf{y}}{\Vert \mathbf{y}-\frac{\mathbf{d}{\left(\mathbf{x}\right)}^{\mathrm{H}}\mathbf{y}\mathbf{d}\left(\mathbf{x}\right)}{{\Vert \mathbf{d}\left(\mathbf{x}\right)\Vert }^2}\Vert }^2,$$

(36)

where y is the actual measurement vector. H represents the conjugate transpose. When the pressure is used as the measurement vector, that is, y = [p(x, ω₁), …, p(x, ω_m), …, p(x, ω_M)]^T, and p(x, ω_m) is the pressure at angular frequency ω_m, with $m\in \left[1,\mathrm{M}\right]$. Furthermore, d(x) = [d₁(x), …, d_m(x), …, d_M(x)]^T, with ${\mathbf{d}}_m\left(\mathbf{x}\right)=2\pi iS\left({\omega }_m\right){\rho }_1{\omega }_m^2{\displaystyle \sum _n{F}_{nm}^2\sin {\beta }_{1nm}H\sin {\beta }_{1nm}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{nm}r\right)}$. When the vertical wave impedance is used as the measurement vector, that is, y = [Z_z(x, ω₁), …, Z_z (x, ω_m), …, Z_z (x, ω_M)]^T, and Z_z (x, ω_m) is the vertical wave impedance. Now, d(x) = [d₁(x), …, d_m(x), …, d_M(x)]^T, with ${\mathbf{d}}_m\left(\mathbf{x}\right)=i{\rho }_1{\omega }_m\frac{{\displaystyle \sum _n{F}_{nm}^2\sin {\beta }_{1nm}H\sin {\beta }_{1nm}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{nm}r\right)}}{{\displaystyle \sum _n{F}_{nm}^2{\beta }_{1nm}\cos {\beta }_{1nm}H\sin {\beta }_{1nm}{z}_s{\mathrm{H}}_0^{\left(1\right)}\left({\xi }_{nm}r\right)}}$. When the measured vectors were the pressure and vertical wave impedance in the frequency domain of 75 Hz–95 Hz, the objective functions for the state vector x are given in Figure 10. When the measured vector was the vertical wave impedance, the range and source depth were very sensitive to the objective function, and the objective functions reached the maximum at 10 km and 40 m. When the measured vector was the pressure, the source depth was sensitive to the objective function. This is why depth estimation could be performed using pressure. However, the range was not sensitive to the objective function, whereby the objective function curve is represented by the black line in Figure 10a. The curve peaked at multiple distances, such as 6 km, 10 km, etc. This multi-peak phenomenon led to instability of ranging performance, as shown in Figure 8a and Figure 9a.

Above all, the source depth will affect the correctness of localization when pressure is used as the measurement vector, while the use of vertical wave impedance to locate the source is not affected by the source depth. In addition, there is no need to know the source spectrum information when using vertical wave impedance for positioning, which is an advantage for passive ranging. In Section 5, the vertical wave impedance is used to locate the explosive sources during an experiment based on the PF.

5. Experimental Results

In 2018, an underwater explosion experiment was conducted in a shallow water waveguide near Qingdao, China. The experimental ship sailed along the scheduled route, and the explosives were placed at a fixed point during the voyage. Two vertical arrays were used as receiving devices. Each vertical array consisted of two hydrophones and one OBS. The navigation trajectory of the experimental ship and the coordinates of the receiving array were obtained from the global positioning system (GPS) data, as shown in Figure 11.

During the experiment, the temperature was recorded by a temperature and depth logger (TD) with a tiny change. Due to the little temperature variability and the small depth of the water column, water sound speed was assumed to be constant. The bottom of the sea could be assumed to be flat, which was estimated as a semi-infinite uniform elastic seabed.

During the experiment, a total of 19 explosives were detonated at four burst spots, and the third and fourth burst spots were selected for sound source localization. Seven explosives (No. 9–No. 15) were selected for source localization. For the seven explosives, the explosive quantity was 50 g. At the third burst spot, four explosives (No. 9–No. 12) were carried out with a depth of 13 m. At the fourth burst spot, three explosives (No. 13–No. 15) were carried out with a depth of 17 m.

Because the second vertical array was on the slope and did not conform to the environment model in this paper, only the data received by the OBS of the first array were used for positioning.

The PF positioning results of experimental data received by the OBS for No. 9–No. 15 explosions are illustrated in Figure 12. The black solid line denotes the estimated value, and the red dashed line denotes the actual distance obtained by the GPS data, with the actual distances shown in

Table 7. Distances of the explosives from the first array and the MAPEs of the localization results.

No.	9	10	11	12	13	14	15
Distance (m)	22,801	22,822	22,844	22,882	30,176	30,158	30,122
MAPE (%) for range	1.11	1.42	1.14	1.14	0.11	0.17	0.29
MAPE (%) for source depth	2.55	2.13	2.86	1.92	3.10	3.30	3.21

. Moreover, Table 7 also gives the MAPE of localization results for No. 9–No. 15 explosions. In the process of using the PF for experimental data, the state process noise covariance value was taken according to the empirical rule obtained in the Section 4. Here, the initial state value was [21000, 15]^T and the ${\mathit{Q}}^{1/2}$ was $\left[\begin{array}{cc}{10}^2& 0\\ 0& 2\end{array}\right]$ for No. 9–No. 12 explosions. Additionally, the initial state value was [29000, 15]^T and the ${\mathit{Q}}^{1/2}$ was $\left[\begin{array}{cc}{10}^2& 0\\ 0& 2\end{array}\right]$ for No. 13–No. 15 explosions. Figure 12 indicates that the experimental range estimates and source depth estimates both converged, and that the differences between the convergence values and the true values were small for No. 9–No. 15 explosions. The MAPEs of range estimation and depth estimation for No. 9–No. 12 did not exceed 1.42% and 2.86%, respectively, while the MAPEs of range estimation and depth estimation for No. 13–No. 15 did not exceed 0.29% and 3.30%. The experimental results proved the accuracy and practicability of the method using vertical wave impedance as a measurement vector based on the PF.

6. Conclusions

In this paper, an OBS sensor was used to estimate the position of the broadband sound source in a Pekeris shallow water waveguide with elastic bottom. In the semi-infinite elastic seabed environment, the expression of pressure and vertical velocity channels received by the OBS sensor were theoretically derived. Based on the PF method, the positioning performances of the pressure and vertical wave impedance as the measurement vector were simulated and analyzed. In the simulation environment, the results showed that the source depth will affect the ability of pressure to locate. When the measured vector was the pressure and the source was 40 m, although the source depth was sensitive to the objective function, the range was not sensitive to the objective function. Moreover, the objective function curve of range peaked at multiple distances. This multi-peak phenomenon led to instability of ranging performance. In contrast, when the measured vector was the vertical wave impedance and the source was 20 m, the range and source depth estimation results converged and were in good agreement with the true values at different frequency bands. In the case of poor localization performance for pressure, that is, when the source was 40 m, the vertical wave impedance as the measurement vector could also exhibit excellent positioning performance. The PF method with the vertical wave impedance as the measurement vector was not affected by source depth and source spectrum information, making it more tolerant and more robust than that with pressure in positioning.

The PF method with the vertical wave impedance as the measurement vector was employed to process the experimental data in the sea near Qingdao. The experimental results showed that the source parameters (the source depth and the range) were correctly estimated and converged to the true values using an OBS sensor in different ranges and different explosive source depths.