Batch Processing through Particle Swarm Optimization for Target Motion Analysis with Bottom Bounce Underwater Acoustic Signals

Abstract

A target angular information in 3-dimensional space consists of an elevation angle and azimuth angle. Acoustic signals propagating along multiple paths in underwater environments usually have different elevation angles. Target motion analysis (TMA) uses the underwater acoustic signals received by a passive horizontal line array to track an underwater target. The target angle measured by the horizontal line array is, in fact, a conical angle that indicates the direction of the signal arriving at the line array sonar system. Accordingly, bottom bounce paths produce inaccurate target locations if they are interpreted as azimuth angles in the horizontal plane, as is commonly assumed in existing TMA technologies. Therefore, it is necessary to consider the effect of the conical angle on bearings-only TMA (BO-TMA). In this paper, a target conical angle causing angular ambiguity will be simulated using a ray tracing method in an underwater environment. A BO-TMA method using particle swarm optimization (PSO) is proposed for batch processing to solve the angular ambiguity problem.

1. Introduction

Acoustic signals are used to indirectly obtain information about objects located underwater. Most passive sonar systems use multiple hydrophones in an array for enhanced performance. A horizontal line array (HLA), used for detecting the azimuth angle of an underwater target, receives acoustic signals with a high signal to noise ratio from designated directions using a beamforming technique. If the target signal intensity is high enough in a designated direction, the target direction is detected. The estimated target direction is represented as a conical angle that indicates the direction of the incoming signal measured by the HLA. Unfortunately, it is impossible to distinguish between up and down or right and left from the conical angle. This is called the cone of ambiguity [1].

Sequential processing and batch processing algorithms are used to estimate the target’s state, including position and velocity, through bearings-only target motion analysis (BO-TMA). There exist several conventional sequential processing algorithms, including the extended Kalman filter [2], the pseudo-measurement filter [3], and the modified gain extended Kalman filter [4]. In addition to these filters, particle filter approaches [5] and random finite set approaches [6,7,8] have been recently introduced. If sufficient computational performance is achieved, sequential processing is suitable for implementation in real time systems. However, good sequential estimation results require small errors in the initial state estimates, and a batch processing algorithm is used for this purpose. Batch processing delivers stable initial values, even though it is not designed to operate in real time because it requires a batch of stored measurements. Robust target localization performance is expected if both types of algorithms are employed properly [9].

In most of the previous studies on sonar systems [2,3,6,7,8,9], it is assumed that the received signal arrives at the HLA through the horizontal plane when the distance between the observer and target is large, or when the observer and the target are located at equal depths. Therefore, the cone of ambiguity of the HLA is simplified to left/right ambiguity, which can be easily addressed through a ship maneuver. However, eigenray tracing results show that the received signal can arrive at the HLA with a high elevation angle, especially along a bottom bounce path [10]. The studies in [11,12] consider the elevation angles in BO-TMA for different sensor depths between the observer and the target. They treat only direct paths without considering the reflection of the ray from the waveguide boundaries (i.e., sea surface and bottom) or the refraction of the ray from the vertical sound speed profile. However, bottom bounce paths, which are generated from the reflection of acoustic waves at the ocean bottom, can produce inaccurate target bearings [13] that affect BO-TMA results.

The ray tracing method [14] is used to calculate the elevation angle due to the refraction and reflection of sound waves in underwater waveguides. This method describes the path of each ray as sound waves propagating through the underwater waveguide. In particular, it is possible to calculate the eigenray [15], which represents the path of a ray that propagates from the source to the receiver. The elevation angle of the target signal can be simulated through eigenray tracing, and the conical angle can be calculated using the azimuth angle and the elevation angle.

In this paper, a study is based on the published conference paper [16] and it is conducted to confirm the observability of TMA using the conical angle including the elevation angle of the path reflected from the bottom interface for a given scenario. A discrete target dynamic equation is established with the target state vector, and the conical angle measurement is obtained from the relative geometry of the observer and the target using the ray tracing method in Section 2. Section 3 presents a method of converting the conical angle into a bearing line in Cartesian coordinates using knowledge of the ocean environment (i.e., bottom bathymetry and a sound speed profile). Additionally, a BO-TMA using the particle swarm optimization (PSO) algorithm is proposed. In Section 4, simulation results for the BO-TMA are analyzed using ray tracing. Finally, a summary and conclusion are given in Section 5.

2. Problem Formulation 2.1. Dynamic Model

The target state vector at the discrete time instancek,1≤k≤K, is defined as:

_Xs(k)=[_pxs(k), _pys(k), _vxs(k), _vys(k)],

_Us(k)=[_uxs(k), _uys(k)],

where_pxs(k)and_pys(k)are the target locations in Cartesian coordinates. Here, thex-axis indicates East and they-axis indicates North. Additionally,_vxs(k)and_vys(k)are the target velocities for each direction, and_uxs(k)and_uys(k)are the target accelerations. The observer state vector is similarly defined as:

_Xo(k)=[_pxo(k), _pyo(k), _vxo(k), _vyo(k)],

_Uo(k)=[_uxo(k), _uyo(k)],

where the subscriptoindicates the observer. Then, the discrete-time system state equation can be described by:

_Xi(k+1)=F_XiT(k)+G_UiT(k),

where_Xiand_Uiare the state vectors of the target (wheni=s) and the observer (wheni=o), and control input, respectively. The superscriptTdenotes a transpose. The state transition matrixFand input coefficient matrixGare defined, respectively, as:

F=[_I2Δt_I202I2],G=[Δ^t22_I2Δt_I2],

where_I2is the 2-dimensional identity matrix,₀₂is the2×2zero matrix, andΔt is the time interval. For system observability, we assume that the sensor outmaneuvers the target while the target is moving with a constant velocity [17].

The horizontal plane trajectories of the target and the observer located at equal depths of 200 m are shown in Figure 1. The total simulation time is 580 s with a sampling period of 20 s so that the total number of scans is 30. The initial state vector of the target,_Xs(1),is[0 m, 2500 m, 0 m/s, −3 m/s]with zero acceleration over the simulation time. The initial state vector of the observer,_Xo(1),is[2000 m, −7000 m, 2.6 m/s, 1.5 m/s]. To ensure system observability, the course of the observer is changed once from 60 to 0° via lateral acceleration starting at 200 s. The bearing change rate is 0.6° per second. The distance between the observer and the target is decreased from a maximum distance of 9.7 km to a minimum distance of 6.9 km.

2.2. Measurement Model

Conventional TMA assumes that the target information obtained from passive line array sonar consists of only the azimuth angle in the horizontal plane, neglecting bottom bounce signals to avoid conical angle ambiguity. In this case, the azimuth angle measured from the north axis,_φn(k), is expressed as:

_φn(k)=atan2(_pxs(k)−_pxo(k), _pys(k)−_pyo(k)),

whereatan2(x, y)denotes a four-quadrant arctangent function that describes the angle between the position of the target and the north axis (positivey-axis). The azimuth angle from the north axis,_φn(k),is converted to the azimuth angle from the direction of the HLA,_φl(k), by subtracting the heading angle of the HLA,_co(k), at each scan timek:

_φl(k)=_φn(k)−_co(k).

In this paper, BO-TMA along with a ray tracing method is used to achieve accurate estimation of target localization in environments with conical angle ambiguity. The conical angle,θ(k), is expressed as:

θ(k)=co^s−1(cos(_φl(k)) ×cos(μ(k)))+v(k),

whereμ(k)is the elevation angle in the vertical plane, andv(k)is the measurement noise modeled as zero mean Gaussian noise with standard deviation_σm. The sign ofθ(k) is unknown from Equation (9), and the conical angle indicates the magnitude of the angle measured from the heading direction of the line array. Thus, the inability to know the exact direction of the arriving signal is known as left/right ambiguity. Various angles used in this paper are shown in Figure 2._φland_φnare the azimuth angles from due north and the direction of the HLA, respectively._co,θ, andμare heading angle of the HLA, conical angle, and elevation angle in the vertical plane, respectively.

A ray tracing method is used to estimate the elevation angle of the target signal in Equation (9). In the ocean, propagation paths of acoustic rays are strongly affected by sound speed profile and bottom bathymetry. These environmental data can be obtained through measurements, from a database, or from an ocean prediction model. In this study, a scenario is constructed that assumes a simple environment. The bathymetry is assumed to be flat with a depth of 2000 m. The sound speed profileC(z) in water is assumed to follow Munk’s sound speed profile and is given by [18]:

C(z)=_C0[1.0+ϵ{^e−η−(1−η)}],

wherezis depth, and_C0is a reference sound speed equal to 1500 m/s as the sound speed at the depth of channel axis_zC(_zC=400 m),η=2(z−_zC)/_zCis a dimensionless depth relative to the channel axis, and the perturbation coefficientϵis equal to7.4×¹⁰⁻³.

The ray paths predicted by the ray tracing method using Munk’s sound speed profile are shown in Figure 3. Although ray tracing was conducted based on observer position, the ray tracing results obtained for opposite directions are the same due to the reciprocity of ray diagrams [14]. In addition, the ray tracing results for all azimuth angles are the same because it is assumed that acoustical ocean parameters are independent of azimuth angle. It is shown in this scenario that only bottom reflected paths exist between the target and the observer, and a direct path from the target does not exist. The elevation angle of the bottom bounce path was calculated by ray tracing to be between 23 and 29° at a target distance of 9.7―6.9 km.

Figure 4 shows the simulation results of the bearing measurements from 30 scans over same time period, which is known as BTR (Bearing-Time Record). The red dashed line that represents the azimuth angle from the north axis_φn(k)was plotted as additional information for assessing the bearing error compared to the conical angle of the bottom bounce path. The bearing error is defined as the difference between_φl(k)and+θ(k)or its mirror angle−θ(k) due to conical angle ambiguity. Figure 4 contains the time histories of_co(k)(the observer heading angle),_φn(k)(the true target azimuth angle),_co(k)+|θ(k)|, and_co(k)−|θ(k)| (two possible bearing angles for TMA that stem from the bottom bounce path). The right/left ambiguity in the horizontal plane is shown in Figure 4 and can be resolved by comparing the histories of_co(k)+|θ(k)|and_co(k)−|θ(k)|. The history of_co(k)−|θ(k)|has smaller variations than that of_co(k)+|θ(k)|. Note that these two angle histories correspond to the history of the true azimuth angle_φn(k). The history of_φn(k) in Figure 4 shows small variations for the entire period that includes the times before and after the observer maneuver, which implies that_co(k)−|θ(k)|rather than_co(k)+|θ(k)|should be applied as the bearing history for this scenario in TMA. From the selection process, the correct sign ofθ(k)in Equation (9) for this scenario is negative. However, even after choosing the bearing history with the correct sign ofθ(k),_co(k)−|θ(k)|still contains bearing error when compared to the true azimuth angle history_φn(k) . Figure 4 shows that this error is ~1° before the observer maneuver and ~13° after the maneuver. This discrepancy is due toμ(k), the elevation angle of the bottom bounce path. Conventional TMA methods for target localization cannot avoid localization errors resulting from bearing errors. Therefore, a new TMA method that accounts for the bottom bounce path is needed.

3. Target Motion Analysis with Bottom Bounce Path 3.1. Bearing Lines of Bottom Bounce Path

Bearing error is due to the elevation angleμ(k)of the bottom bounce path, which is unknown even after selection of the correct sign ofθ(k). In this study, the bearing line in Cartesian coordinates is introduced. Define thei-th expected azimuth angle_φ^l(k, i)for1≤i≤I, which represents a possible target azimuth angle relative to the heading direction of the HLA. According to Equation (9),_φ^l(k, i)must lie within the range between zero and the conical angleθ(k), and then the elevation angleμ^(k, i)can be estimated as:

μ^(k, i)=co^s−1(cos(θ(k))cos(_φ^l(k, i))).

The sign of_φ^l(k, i)is equal to the sign ofθ(k). For each_φ^l(k, i), ray tracing for the ray launched at an angle ofμ^(k, i)from the observer position is conducted to find the ranger^(k, i)of target location if it exists in the direction of_φ^l(k, i) (Figure 5a). Since the target depth was assumed to be 200 m, the distance at which the ray arrives at a water depth of 200 m after bottom reflection becomes the target range in the_φ^l(k, i)direction. This process is repeatedi=I times (Figure 5b). In this study, the expected azimuth angle was varied every 0.5°. Accordingly,|θ(k)|divided by 0.5° was used to determine the value ofIfor each scank.

For thek-th scan,Ipossible target positions in the horizontal plane corresponding to every_φ^l(k, i)are connected in a line, which is defined as a bearing line in this paper. The possible target position vector in Cartesian coordinates with_φ^l(k, i)andr^(k, i)is denoted as:

L^(k, i)=[_p^xl(k, i), _p^yl(k, i)].

Figure 6 is drawn in the horizontal plane and it shows the bearing lines corresponding tok=1andk=K. The lines (denoted by line of conical angle) indicating the measured conical angleθ(k)in the horizontal plane fork=1andk=K. If the elevation angle is not considered, as in previous studies, the bearing line is displayed as a straight line. However, the bearing lineL^(k, i)is displayed as a curved line when the elevation angle is considered. Conventional batch estimation methods for TMA utilize the conical angles to determine the initial target states, while the proposed TMA method utilizes the bearing lines. The objective of the proposed TMA method is to find the optimal initial position and velocity of the target based on the bearing lines in Cartesian coordinates using the PSO algorithm to minimize the objective function.

3.2. Particle Swarm Optimization

The PSO algorithm is a stochastic optimization algorithm used to find the optimal positions of particles and is based on the social behavior of animals moving in flocks [19,20]. In BO-TMA studies, each particle representing an estimated initial target state vector consists of four elements: the positions and velocities in thexandydirections. First, atk=1, the particles are uniform, randomly spread along the bearing line within the target-observer distance from 1 to 30 km. A specific velocity vector, which is randomly selected in the range of0|v^| 10m/s, where|v^|=^{|_v^x|2}+^{|_v^y|2}, is assigned to each particle. Then, the position of each particle at the next scan time (k=2) is calculated using the dynamic model shown in Section 2.1 from the position atk=1. In this manner, a total ofKpositions are determined for each particle, which forms a particle trajectory. After that, the shortest distance between each particle position and the bearing line corresponding to the same scan time number is calculated. This distance is then normalized by the distance between the observer and particle position at each scan time to avoid excessive convergence to local optima, which happens because distance error increases as the distance between the observer and the particle increases. Finally, the normalized distance errors for allKparticle positions are summed to obtain an objective function^Jmfor them-th particle, which is expressed as:

^Jm=∑k=1Kmini^{(_p^xm(k)−_p^xl(k, i))2}+^{(_p^ym(k)−_p^yl(k, i))2(_p^xm(k)−_pxo(k))2}+^{(_p^ym(k)−_pyo(k))2},

where_p^xm(k)and_p^ym(k), respectively, are the positions in thexandydirections of them-th particle at scan timek; and_pxo(k)and_pyo(k)are the observer positions in thexandydirections at scan timek . The total particle number used here was 200 (Table 1). Since each particle is considered a candidate for the target, the next step is to find the initial state vector of the particle that produces the minimum value of^Jm . In this study, the PSO algorithm was used as an optimization technique to find the optimal target trajectory. In each generation, the best values for the state vectors consisting of the positions and velocities of the particles are evaluated by comparison with state vectors selected during previous generations, and the state change rates of the particles are adjusted based on the experiences of the particles and their companions. The state vectors in the next generation are updated with the sum of the present state vectors and the adjusted state change rates of the particles [20]. The process is expressed as [19]:

_vp(n+1, m, d)=_{c1 vp}(n,m,d)+_vl(n,m,d)+_vs(n,m,d),

_vl(n,m,d)=_{c2 r1}{_xpl(m,d)−_xp(n,m,d)},

_vs(n,m,d)=_{c3 r2}{_xps(d)−_xp(n,m,d)},

and _xp(n+1,m,d)=_xp(n,m,d)+_vp(n+1, m, d),

where_xp(n,m, d)represents the state vector of them-th particle for then-th generation with dimensiond. Dimensiondis one of 1, 2, 3, and 4 corresponding to the positions and velocities of the particles atk = 1, that is,_p^xm(1),_p^ym(1),_v^xm(1), and_v^ym(1), respectively. In addition,_vp(n,m, d)represents the state change rates of the particles for_xp(n,m, d). Finally,_vl(n,m, d)and_vs(n,m, d)are the local state change rate and the social state change rate for them-th particle, respectively. The local state vector_xpl(m, d)is the best state vector of them-th particle obtained from the first generation to then-th generation, and the social state vector_xps(d)is the best state vector of the particle with the smallest^Jmof all particles up to then-th generation. In the above equations,_c1,_c2, and_c3are acceleration weight constants determined empirically through many trial runs to be 0.73, 0.1, and 0.2, respectively. Random numbers_r1and_r2are selected in the range between 0 and 1. The process is iterated until the state vector of each particle converges to the best state vector that satisfies the minimum position errors. In our case, the generation is terminated when the standard deviations of the positions,_σp, and velocities,_σv, of 200 particles converge to values less than 100 m and 0.2 m/s, respectively. Finally, the trajectory of the particle with the best state vector is selected as the target trajectory.

4. Simulation Result

For the observability test, it was assumed that the water depth was 2000 m and the bottom topography was flat. The conical angle was calculated using the azimuth and elevation angle of the acoustic ray path between the target and the observer. Munk’s sound speed profile was used for ray tracing to calculate the elevation angle. To test the applicability of batch processing using the PSO algorithm proposed in this paper, it was assumed that Gaussian noise with zero mean and standard deviation_σmwas included in the conical angle measurements. Three values of_σm (0.2, 0.4, and 0.6°) were considered for comparison purposes. For this scenario, the conical angle was estimated to change at a rate of approximately 0.5°/scan except during the period in which the observer heading changed. Figure 7 shows the histories of conical angles with measurement errors corresponding to three different standard deviations.

One thousand random runs were generated for each of the three standard deviations of the conical angle measurement errors, and TMA was carried out for each run. The results are shown in Figure 8, in which the left column shows the scatter plots for the estimates of initial target position for the 1000 runs, and the right column shows the scatter plots of target velocity. For the different standard deviations, the mean values of the estimated initial state vectors and their variances are listed in Table 2. The results show that, as the standard deviation of the measurement error increases, the distribution of the initial state vector obtained from the proposed BO-TMA becomes wider. In particular, as the measurement error increases, the estimated positions of the target tend to spread wider along the bearing line atk=1, which is reasonable because the particles were spread along the bearing line atk=1. The mean value of the initial state vector estimated for the standard deviation of 0.2° (marked by a yellow triangle in the figure) has the best agreement with the true initial target state vector (marked by red circle), and as the standard deviation increases, the difference increases slightly. However, the mean values for the three cases are still in good agreement with the true values.

To investigate the accuracy of the TMA results with increasing the number of scansk , the processes were repeated with the scan numbers of 15, 30, and 60 which correspond to the sampling periods of 40, 20, and 10 s, respectively. The standard deviation of the conical angle measurements were assumed to be 0.4°. The estimation results of the initial target state vector with the three scan numbers are shown in Figure 9, and the resulting mean values and variances are listed in Table 3. Figure 9 and Table 3 indicate that more frequent collection of conical angle measurements achieves more accurate TMA results with increased expense of computational resources.

As shown in Figure 7, the bearing errors due to elevation angle after the observer maneuver are larger than 10°. Conventional TMA methods produce large localization errors in environments dominated by acoustic rays being strongly reflected or refracted up and down. However, the proposed BO-TMA method using ray tracing shows good localization performance in such environments, which implies that the proposed TMA method is a more effective tool for increasing solution accuracy in real underwater applications, especially in waveguide environments where bottom bounce paths are dominant.

5. Summary and Conclusion In this paper, a BO-TMA algorithm using a ray tracing method is proposed to accurately consider the conical angles generated by bottom bounce paths. The 3-dimensional conical angle information was converted to bearing lines in a 2-dimensional plane using a ray tracing method. Then, the PSO algorithm was carried out based on the constructed bearing lines to find optimal target state vectors. The BO-TMA method using ray tracing and the PSO algorithm proposed in this paper is summarized below.

(1) Convert the conical angles of the bottom bounce path to a bearing line using the ray tracing technique. Set the generation numbern=1.

(2) Initialize particles with the bearing line atk=1. Uniform, randomly spread particles on the bearing line and assign velocities randomly selected in the range0≤|v^|≤10m/s.

(3) For each particle with a four-element state vector, calculate the objective function^Jmusing the particle trajectories and the bearing lines corresponding tok=1,⋯,K.

(4) Find the particle that produces the minimum value of^Jm.

(5) Generate the next generation particle group by applying the PSO algorithm. (6) Go to Step (3), and then iterate the process. (7) Terminate the iteration when the state vectors of the particles reach the termination condition. In this paper, a ray tracing technique was used to calculate the elevation angle. The conical angle of the target was then calculated based on the estimated elevation angle. Characteristics of the oceanic environment are known, allowing for accurate estimation of elevation angles. However, since the oceanic environment fluctuates temporally and spatially, errors can arise from uncertainty in environmental information. In addition, we assumed that target depth is the same as observer depth. Uncertainty in target depth may also result in target distance errors. Therefore, further research into various target-observer geometries and various ocean environments is required to generalize the results shown in this paper.

Figure 1. Trajectories of (a) the target and (b) the observer in the horizontal plane.

Figure 2. Geometry between observer and target.φlandφnare the azimuth angles from due north and the direction of the HLA, respectively.co,θ, andμare heading angle of the HLA, conical angle, and elevation angle in the vertical plane, respectively.

Figure 3.(a) Ray paths predicted by ray tracing method based on (b) Munk's sound speed profile. Direct and bottom bounce paths are plotted with magenta and blue lines, respectively.

Figure 4. Bearing-time records (BTRs) of the scenario.

Figure 5.(a) Eigenray tracing result conducted to determine the expected target range. The distance at which the ray arrives at an expected target depth after bottom reflection becomes the estimated target range inφ^l(k, i) direction. (b) Top-view illustration showing the line of conical angle and bearing line. For k-th scan, the line connecting I possible target positions estimated using the eigenray tracing is a bearing line (red line in figure).

Figure 6. Bearing lines (solid lines) and lines of conical angles (dashed lines) atk=1andk=K.

Figure 7. The BTRs for conical angle measurements including Gaussian measurement error with zero mean and standard deviation of (a) 0.2, (b) 0.4, and (c) 0.6°.

Figure 8. The distribution of initial states estimated using TMA for 1000 random runs for standard deviations of zero mean Gaussian measurement errors of (a) 0.2, (b) 0.4, and (c) 0.6°. The true initial state vector of the target is[0 m, 2500 m, 0 m/s, -3 m/s]. The left column shows the initial target position estimates, and the right shows target velocity estimates. The true initial state vector of the target and the mean of estimated state vectors are indicated by red circles and yellow triangles, respectively. The regions within one standard deviation of the mean are indicated by black ellipses.

Figure 8. The distribution of initial states estimated using TMA for 1000 random runs for standard deviations of zero mean Gaussian measurement errors of (a) 0.2, (b) 0.4, and (c) 0.6°. The true initial state vector of the target is[0 m, 2500 m, 0 m/s, -3 m/s]. The left column shows the initial target position estimates, and the right shows target velocity estimates. The true initial state vector of the target and the mean of estimated state vectors are indicated by red circles and yellow triangles, respectively. The regions within one standard deviation of the mean are indicated by black ellipses.

Figure 9. The distribution of initial states estimated using TMA for 1000 random runs with standard deviations of zero mean Gaussian measurement error of 0.4° with the measurement numbers of (a) 15, (b) 30, and (c) 60.

Parameter	Symbol	Value
Number of particles	m	200
Number of dimensions	d	4
Number of generations	n	σp<100m andσv<0.2m/s
Acceleration weight constantsc1	c1	0.73
Acceleration weight constantsc2	c2	0.1
Acceleration weight constantsc3	c3	0.2
Random numberr1	r1	0―1
Random numberr1	r2	0―1

Standard Deviation of Measurement Noise, σm	Mean of Initial Target State Vector, Xs^,[m, m, m/s, m/s]	Variance of Initial Target State Vector, σs^2,[m2, m2, m/s2, m/s2]
0.2°	[18, 2524, −0.2, −2.8]	[712, 1782, 0.52, 0.52]
0.4°	[5, 2607, −0.4, −2.9]	[1052, 3402, 0.62, 0.72]
0.6°	[−6, 2693, −0.5, −2.8]	[1382, 4732, 0.72, 0.92]

Number of Measurements, k	Mean of Initial Target State Vector, Xs^,[m, m, m/s, m/s]	Variance of Initial Target State Vector, σs^2,[m2, m2, m/s2, m/s2]
15	[−18, 2655, −0.2, −3.0]	[1382, 5092, 0.62, 0.92]
30	[5, 2607, −0.4, −2.9]	[1052, 3402, 0.62, 0.72]
60	[11, 2590, −0.4, −2.9]	[822, 2362, 0.52, 0.62]

Author Contributions

Conceptualization, R.O. and J.W.C.; methodology, R.O. and T.L.S.; software, R.O.; validation, R.O., T.L.S. and J.W.C.; formal analysis, R.O. and T.L.S.; investigation, R.O.; resources, R.O.; data curation, R.O.; writing-original draft preparation, R.O.; writing-review and editing, T.L.S. and J.W.C.; visualization, R.O.; supervision, T.L.S. and J.W.C.; project administration, J.W.C.; funding acquisition, J.W.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was a part of the project titled "Development of Distributed Underwater Monitoring and Control Networks", funded by the Ministry of Oceans and Fisheries, Korea.

Conflicts of Interest

The authors declare no conflict of interest.

References

1. Abraham, D.A. Underwater Acoustic Signal Processing: Modeling, Detection, and Estimation, Modern Acoustics and Signal Processing; Springer Nature Switzerland AG: Cham, Switzerland, 2019.

2. Aidala, V.J.; Hammel, S.E. Utilization of Modified Polar Coordinates for Bearings-Only Tracking. IEEE Trans. Autom. Control 1983, AC-28, 283-294.

3. Aidala, V.J.; Nardone, S.C. Biased Estimation Properties of the Pseudolinear Tracking Filter. IEEE Trans. Aerosp. Electron. Syst. 1982, AES-18, 432-441.

4. Song, T.L.; Speyer, J.L. A stochastic analysis of a modified gain extended Kalman filter with applications to estimation with bearings only measurements. IEEE Trans. Autom. Control 1985, AC-30, 940-949.

5. Arulampalam, S.; Maskell, S.; Gordon, N.; Clapp, T. A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Trans. Signal process. 2002, 50, 174-188.

6. Beard, M.; Vo, B.T.; Vo, B.N.; Arulampalam, S. A partially uniform target birth model for gaussian mixture PHD/CPHD Filtering. IEEE Trans. Aerosp. Electron. Syst. 2013, 49, 2835-2844.

7. Zhang, Q.; Song, T.L. Improved bearings-only multi-target tracking with GM-PHD filtering. Sensors 2016, 16, 1469.

8. Xie, Y.; Song, T.L. Bearings-only multi-target tracking using an improved labeled multi-Bernoulli filter. Signal Process. 2018, 151, 32-44.

9. Kirubarajan, T.; Shalom, Y.B.; Lerro, D. Bearings-Only Tracking of Maneuvering Targets Using a Batch-Recursive Estimator. IEEE Trans. Aerosp. Electron. Syst. 2001, 37, 770-780.

10. Oh, R.; Gu, B.S.; Kim, S.; Song, T.L.; Choi, J.W.; Son, S.U.; Kim, W.K.; Bae, H.S. Estimation of bearing error of line array sonar system caused by bottom bounced path. J. Acoust. Soc. Korea 2018, 37, 412-421.

11. Jauffret, C.; Pillon, D. Observability in Passive Target Motion Analysis. IEEE Trans. Aerosp. Electron. Syst. 1996, 32, 1290-1300.

12. Hammel, S.E.; Aidala, V.J. Observability Requirements for Three-Dimensional Tracking via Angle Measurements. IEEE Trans. Aerosp. Electron. Syst. 1985, AES-21, 200-207.

13. Heaney, K.D.; Campbell, R.L.; Murray, J.J. Deep water towed array measurements at close range. J. Acoust. Soc. Am. 2013, 134, 3230-3241.

14. Jensen, F.B.; Kuperman, W.A.; Porter, M.B.; Schmidt, H. Computational Ocean Acoustics, AIP Series in Modern Acoustics and Signal Processing; American Institute of Physics: New York, NY, USA, 1994.

15. Porter, M.B.; Bucker, H.P. Gaussian beam tracing for computing ocean acoustic fields. J. Acoust. Soc. Am. 1987, 82, 1349-1359.

16. Oh, R.; Song, T.L.; Choi, J.W. Bearings-Only Target Motion Analysis using Ray Tracing. In Proceedings of the 8th International Conference on Control Automation & Information Sciences, Chengdu, China, 23-26 October 2019.

17. Song, T.L. Observability of target tracking with bearings-only measurements. IEEE Trans. Aerosp. Electron. Syst. 1996, 32, 1468-1472.

18. Munk, W.H. Sound channel in an exponentially stratified ocean with application to SOFAR. J. Acoust. Soc. Am. 1974, 55, 220-226.

19. Kennedy, J. Particle Swarm Optimization. In Encyclopedia of Machine Learning; Springer: New York, NY, USA, 2011; pp. 760-766.

20. Kim, S.; Choi, J.W. Optimal Deployment of Sensor Nodes Based on Performance Surface of Underwater Acoustic Communication. Sensors 2017, 17, 2389-2403.