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

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.


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 Sensors 2020, 20, 1234 2 of 14 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.

Dynamic Model
The target state vector at the discrete time instance k, 1 ≤ k ≤ K, is defined as: where p xs (k) and p ys (k) are the target locations in Cartesian coordinates. Here, the x-axis indicates East and the y-axis indicates North. Additionally, v xs (k) and v ys (k) are the target velocities for each direction, and u xs (k) and u ys (k) are the target accelerations. The observer state vector is similarly defined as: where the subscript o indicates the observer. Then, the discrete-time system state equation can be described by: where X i and U i are the state vectors of the target (when i = s) and the observer (when i = o), and control input, respectively. The superscript T denotes a transpose. The state transition matrix F and input coefficient matrix G are defined, respectively, as: where I 2 is the 2-dimensional identity matrix, 0 2 is the 2 × 2 zero 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, X s (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, X o (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. where and are the state vectors of the target (when = ) and the observer (when = ), and control input, respectively. The superscript denotes a transpose. The state transition matrix and input coefficient matrix are defined, respectively, as: where 2 is the 2-dimensional identity matrix, 0 2 is the 2 × 2 zero matrix, and ∆ 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, (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, (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.

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, ( ), is expressed as: where 2( , ) denotes a four-quadrant arctangent function that describes the angle between the position of the target and the north axis (positive -axis). The azimuth angle from the north axis, ( ), is converted to the azimuth angle from the direction of the HLA, ( ), by subtracting the heading angle of the HLA, ( ), at each scan time : 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, ( ), is expressed as:

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: where atan2(x, y) denotes a four-quadrant arctangent function that describes the angle between the position of the target and the north axis (positive y-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, c o (k), at each scan time 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) = cos −1 (cos(ϕ l (k)) × cos(µ(k))) + v(k), where µ(k) is the elevation angle in the vertical plane, and v(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. ϕ l and ϕ n are the azimuth angles from due north and the direction of the HLA, respectively. c o , θ, and µ are heading angle of the HLA, conical angle, and elevation angle in the vertical plane, respectively. where ( ) is the elevation angle in the vertical plane, and ( ) is the measurement noise modeled as zero mean Gaussian noise with standard deviation . The sign of ( ) 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. and are the azimuth angles from due north and the direction of the HLA, respectively. , , and are heading angle of the HLA, conical angle, and elevation angle in the vertical plane, respectively. and are the azimuth angles from due north and the direction of the HLA, respectively. , , 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 profile ( ) in water is assumed to follow Munk's sound speed profile and is given by [18]: where is depth, and 0 is a reference sound speed equal to 1500 m/s as the sound speed at the depth of channel axis ( = 400 m), = 2( − )/ is a dimensionless depth relative to the channel axis, and the perturbation coefficient is equal to 7.4 × 10 −3 .
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. 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 profile C(z) in water is assumed to follow Munk's sound speed profile and is given by [18]: where z is depth, and C 0 is a reference sound speed equal to 1500 m/s as the sound speed at the depth of channel axis z C (z C = 400 m), η = 2(z − z C )/z C is a dimensionless depth relative to the channel axis, and the perturbation coefficient is equal to 7.4 × 10 −3 . 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 ( ) 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 ( ) and + ( ) or its mirror angle − ( ) due to conical angle ambiguity. Figure 4 contains the time histories of ( ) (the observer heading angle), ( ) (the true target azimuth angle), ( ) + | ( )|, and ( ) − | ( )| (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 ( ) + | ( )| and ( ) − | ( )|. The history of ( ) − | ( )| has smaller variations than that of ( ) + | ( )|. Note that these two angle histories correspond to the history of the true azimuth angle ( ). The history of ( ) in Figure  4 shows small variations for the entire period that includes the times before and after the observer maneuver, which implies that ( ) − | ( )| rather than ( ) + | ( )| should be applied as the bearing history for this scenario in TMA. From the selection process, the correct sign of ( ) in Equation (9) for this scenario is negative. However, even after choosing the bearing history with the correct sign of ( ), ( ) − | ( )| still contains bearing error when compared to the true azimuth angle history ( ). Figure 4 shows that this error is ~1° before the observer maneuver and ~13° after the maneuver. This discrepancy is due to ( ), 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.  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 c o (k) (the observer heading angle), ϕ n (k) (the true target azimuth angle), c o (k) + θ(k) , and c o (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 c o (k) + θ(k) and c o (k) − θ(k) . The history of c o (k) − θ(k) has smaller variations than that of c o (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 c o (k) − θ(k) rather than c o (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), c o (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.

Bearing Lines of Bottom Bounce Path
Bearing error is due to the elevation angle ( ) of the bottom bounce path, which is unknown even after selection of the correct sign of ( ). In this study, the bearing line in Cartesian coordinates is introduced. Define the -th expected azimuth angle � ( , ) for 1 ≤ ≤ , which represents a possible target azimuth angle relative to the heading direction of the HLA. According to Equation (9), � ( , ) must lie within the range between zero and the conical angle ( ), and then the elevation angle ( , ) can be estimated as: The sign of � ( , ) is equal to the sign of ( ). For each � ( , ), ray tracing for the ray launched at an angle of ( , ) from the observer position is conducted to find the range ( , ) of target location if it exists in the direction of � ( , ) (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 � ( , ) direction. This process is repeated = times (Figure 5b). In this study, the expected azimuth angle was varied every 0.5°. Accordingly, | ( )| divided by 0.5° was used to determine the value of for each scan .

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 the i-th expected azimuth angleφ l (k, i) for 1 ≤ 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) = cos −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 repeated i = 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 of I for each scan k.
For the k-th scan, I possible 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: Figure 6 is drawn in the horizontal plane and it shows the bearing lines corresponding to k = 1 and k = K. The lines (denoted by line of conical angle) indicating the measured conical angle θ(k) in the horizontal plane for k = 1 and k = 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  If the elevation angle is not considered, as in previous studies, the bearing line is displayed as a straight line. However, the bearing line � ( , ) 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.    If the elevation angle is not considered, as in previous studies, the bearing line is displayed as a straight line. However, the bearing line � ( , ) 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.

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 the x and y directions. First, at k = 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 Sensors 2020, 20, 1234 8 of 14 is randomly selected in the range of 0|v| 10 m/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 at k = 1. In this manner, a total of K positions 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 all K particle positions are summed to obtain an objective function J m for the m-th particle, which is expressed as: wherep xm (k) andp ym (k), respectively, are the positions in the x and y directions of the m-th particle at scan time k; and p xo (k) and p yo (k) are the observer positions in the x and y directions at scan time k.
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 J m . 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]: and x p (n + 1, m, d) = x p (n, m, d) + v p (n + 1, m, d), where x p (n, m, d) represents the state vector of the m-th particle for the n-th generation with dimension d. Dimension d is one of 1, 2, 3, and 4 corresponding to the positions and velocities of the particles at k = 1, that is,p xm (1),p ym (1),v xm (1), andv ym (1), respectively. In addition, v p (n, m, d) represents the state change rates of the particles for x p (n, m, d). Finally, v l (n, m, d) and v s (n, m, d) are the local state change rate and the social state change rate for the m-th particle, respectively. The local state vector x pl (m, d) is the best state vector of the m-th particle obtained from the first generation to the n-th generation, and the social state vector x ps (d) is the best state vector of the particle with the smallest J m of all particles up to the n-th generation. In the above equations, c 1 , c 2 , and c 3 are acceleration weight constants determined empirically through many trial runs to be 0.73, 0.1, and 0.2, respectively. Random numbers r 1 and r 2 are 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. Table 1. Particle swarm optimization parameters used to find the optimal initial position and velocity of target.

Parameter Symbol Value
Number

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 σ m was 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.
Sensors 2020, 20, 1234 9 of 14 deviations of the positions, , and velocities, , 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. Table 1. Particle swarm optimization parameters used to find the optimal initial position and velocity of target.

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 was included in the conical angle measurements. Three values of (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 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 at k = 1, which is reasonable because the particles were spread along the bearing line at k = 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.  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 at = 1, which is reasonable because the particles were spread along the bearing line at = 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 scans , 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. To investigate the accuracy of the TMA results with increasing the number of scans k, 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.

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 number n = 1.
(2) Initialize particles with the bearing line at k = 1. Uniform, randomly spread particles on the bearing line and assign velocities randomly selected in the range 0 ≤ |v| ≤ 10 m/s.
(3) For each particle with a four-element state vector, calculate the objective function J m using the particle trajectories and the bearing lines corresponding to k = 1, · · · , K.
(4) Find the particle that produces the minimum value of J m . (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.