Three-Dimensional Tracking of a Target under Angle-Frequency Measurements with Multiple Frequency Lines

This article considers tracking a constant-velocity underwater target, which emits sound with distinct frequency lines. By analyzing the target’s azimuth, elevation and multiple frequency lines, the ownship can estimate the target’s position and (constant) velocity. In our paper, this tracking problem is called the 3D Angle-Frequency Target Motion Analysis (AFTMA) problem. We consider the case where some frequency lines disappear and appear occasionally. Instead of tracking every frequency line, this paper proposes to estimate the average emitting frequency by setting the average frequency as the state vector in the filter. As the frequency measurements are averaged, the measurement noise decreases. In the case where we use the average frequency line as our filter state, both the computational load and the root mean square error (RMSE) decrease, compared to the case where we track every frequency line one by one. As far as we know, our manuscript is unique in addressing 3D AFTMA problems, such that an ownship can track an underwater target while measuring the target’s sound with multiple frequency lines. The performance of the proposed 3D AFTMA filter is demonstrated utilizing MATLAB simulations.


Introduction
This article considers a scenario where an ownship tracks a constant-velocity underwater target, which emits sound with distinct frequency lines. For instance, this scenario can be utilized for tracking of an enemy submarine by processing the sound generated from the submarine.
Multi-frequency lines, concurrently processed with bearing measurements, can be used for tracking a constant-velocity target [1]. For example, the target's engines or propellers make noise with distinct frequency lines. Moreover, some frequency lines may disappear and appear occasionally [2].
Inspired by [1], our paper handles the case where a constant-velocity target generates multiple frequency lines. Reference [1] addressed a typical architecture for processing underwater sonar measurements, so that a target's bearing and frequency lines can be estimated in real time.
By analyzing the target's azimuth, elevation and multiple frequency lines, the ownship can estimate the target's position and (constant) velocity. In our paper, this 3D tracking problem is called the 3D Angle-Frequency Target Motion Analysis (AFTMA) problem. The proposed tracking filter can handle multiple frequency lines as well as disappearance of a frequency line. As far as we know, no paper in the literature on 3D AFTMA has considered both multiple frequency lines and disappearance of a frequency line.
Considering 2D scenarios, the ownship can track a target based on bearing and frequency measurements [1,[3][4][5][6][7]. Considering 2D cluttered environments, ref. [4] addressed a solution to multi-target tracking based on the Gaussian mixture probability hypothesis density (GM-PHD) filter jointly using the frequency and bearing measurements. In the case where the ownship can measure the target's frequency and bearing, the ownship's maneuver is not necessary for tracking a constant-velocity target, provided that the unknown emission frequency is a constant and that the bearing rate is not equal to zero [1].
Considering 3D scenarios, refs. [8,9] addressed the 3D AFTMA problem, where the ownship tracks a target based on frequency, elevation, and bearing measurements. The authors of [9] addressed the 3D AFTMA, by considering measurement parameters such as frequency, elevation and bearing of the target. In ref. [8], the Gauss-Helmert model was introduced to solve the 3D AFTMA problem by expanding the unknown signal emission time as an unknown variable to the state. The papers reviewed in this paragraph did not consider the case where a target generates multiple frequency lines. Note that this distinguishes our paper from other papers on 3D AFTMA.
Reference [1] considered 2D scenarios for tracking a constant-velocity target generating multiple frequency lines. However, ref. [1] did not consider 3D environments, as in our paper. In addition, ref. [1] did not consider the fact that a frequency line may disappear and appear occasionally. Ref. [1] considered the ideal case where all frequency lines exist at all times. However, ref. [1] cannot handle a practical scenario where a frequency line may disappear and appear occasionally.
Reference [2] considered 2D scenarios for tracking a constant-velocity target, in the case where there are two frequency lines. Reference [2] considered the case where a frequency line appears and disappears occasionally. However, ref. [2] considered 2D scenarios where only two frequency lines exist. Distinct from [1,2], our paper addresses the 3D AFTMA problem, considering the case where multiple frequency lines can disappear and appear occasionally.
Instead of tracking every frequency line one by one, we propose to estimate the average emitting frequency by setting the average frequency as the state vector in the filter. As the frequency measurements are averaged, the measurement noise decreases. Moreover, using the average frequency as our Kalman filter state, we can decrease the matrix size of the covariance matrix in the Kalman filter. This decreases the numerical error in matrix calculations, such as the calculation of an inverse matrix. Using MATLAB simulations, we show that as we use the average frequency line as our filter state, both the computational load and the root mean square error (RMSE) decrease, compared to the case where we track every frequency line one by one.
As far as we know, our paper is novel in addressing the 3D AFTMA with multiple frequency lines, such that a frequency line may disappear or appear occasionally. The performance of the proposed AFTMA filter is further demonstrated utilizing MATLAB simulations.
This paper is organized as follows. Section 2 presents the 3D AFTMA problem handled in this paper. Section 3 presents the proposed tracking filter for solving the 3D AFTMA problem with multiple frequency lines. Section 4 presents the MATLAB simulations. Section 5 presents the conclusions.

3D AFTMA Problem
Before addressing the 3D AFTMA problem, some definitions are presented. Let I a denote an identity matrix with size a × a. Let 0 a denote a zero matrix with size a × a. Let 0 a,1 denote a zero vector with size a × 1. In addition, diag(a, b, c, . . .) denotes a diagonal matrix with diagonal elements a, b, c, . . . in this order. Let C[i] denote the i-th element in a column vector C. Let C[1 : n] denote a column vector composed of first n elements in C. Considering a matrix M, M[1 : n, 1 : n] indicates the sub-matrix of M, consisting of the first n rows and n columns of M. Let Cov(q) denote the error covariance for a variable q for convenience. Let N(0, M) indicate the Gaussian distribution with mean 0 and variance M.
This manuscript addresses discrete time systems. Let dt indicate the sampling interval in discrete time systems. At sample-step k, some frequency lines are measured by the ownship [1]. Suppose that M k frequency lines are measured at sample-step k. At samplestep k, we sort measured frequency lines in ascending order. for all m ∈ {1, 2, . . . , M k }. Since frequency lines can appear and disappear as time goes on, S e k = S e k is feasible for k = k . It is assumed that the emitting frequency of a frequency line is a constant. This assumption has been widely applied in tracking problems based on frequency measurements [1][2][3][4]8,10]. For instance, suppose that an emitting frequency, say f e , exists at all sample-steps k. Then, f e ∈ S e k for all sample-steps k. For all m ∈ {1, 2, . . . , M k }, f e,m k is not known a priori. Suppose that one builds a filter of estimating S e k , i.e., we wish to estimate f e,m k for all m ∈ {1, 2, . . . , M k }. This implies that we build a filter for tracking every frequency line whenever a new line appears.
If M k increases as sample-step k increases, then the number of unknown variables in S e k increases in the filter. This increases the size of the state vector in the filter. Thus, increasing M k increases the computational load for estimation of unknown emitting frequencies in S e k . Instead of estimating all emitting frequencies in S e k , the proposed approach is to estimate the average emitting frequency by setting the average frequency as the state vector in the filter. In other words, we propose to estimate the average emitting frequency by setting the average frequency as the state vector in the filter.
The authors of [1] considered the ideal case where all frequency lines exist at all times. For this ideal case, we derived the following results under MATLAB simulations (see Section 4). As we track every frequency line one by one, the computational load increases, compared to the case where we use the average frequency line as our filter state. Moreover, as we track every frequency line one by one, the RMSE increases, compared to the case where we use the average frequency line as our filter state.
By setting the average frequency as the state vector in the filter, increasing M k does not increase the size of the state vector in the filter. This study thus defines the average emitting frequency as f e,avg k Note that f e,avg k in (3) is not known in advance, since f e,m k is not known a priori. Hence, one addresses a time-efficient filter for real-time estimation of f e,avg k . Recall that the sorted emitting frequency set S e k is defined in (2). In practice, there may be a case where S e k−1 = S e k . This implies that a frequency line can appear or disappear at sample-step k. A frequency line change at sample-step k leads to a sudden jump in the average emitting frequency f e,avg k .

Process Model
The target's 3D position, velocity, and emitting frequency at sample-step k are denoted as Considering a constant velocity (constant speed and heading) target, the target's motion model is where n k = Γ × (a x , a y , a z ) T . In n k , (a x , a y , a z ) T denotes a target's acceleration that is not known in advance. In addition, Γ is In (4), F is In (4), in which σ a presents the standard deviation of a. Our tracking filter uses the relative state vector X k as follows.
) T denote the state vector in our tracking filter. The relative state vector X k is distinct from the target state vector X E k in (4). Under (4) and (8), one derives

Measurement Model
In the 3D AFTMA problem, the target's elevation, azimuth, and frequency are measured at sample-step k. Suppose that the ownship measures M k frequency lines at samplestep k. Let f m k denote the m-th frequency line measurement (m ∈ {1, 2, . . . , M k }). The frequency measurement equation is derived as follows. Let C denote the underwater sound speed, which is assumed to be accessible a priori. This assumption has been widely applied in tracking problems based on frequency measurements [1][2][3][4]8,10]. Then, In (10), n f is the measurement noise with Gaussian distribution given as n f ∼ N(0, σ 2 f ). We assume that the standard deviation of every frequency line is identical.
In (10), r k = X k [1 : 3] denotes the range between the ownship and the target at sample-step k. Thus,ṙ k in (10) iṡ The entire measurements at sample-step k are represented as In (12), we use and Here, b 1 (X k ) and b 2 (X k ) define the azimuth and elevation angles, respectively. In (13), n b 1 denotes the azimuth measurement noise with Gaussian distribution given as n b 1 ∼N(0, σ 2 b,1 ). In (14), n b 2 denotes the elevation measurement noise with Gaussian distribution given as . Using (12) as the measurement equation, we can build a 3D AFTMA filter for tracking every frequency line whenever a new frequency line appears. However, as the number of frequency lines increases, the state vector size of the filter increases. The authors of [1] stated that more than 100 frequency lines can exist in practice. In the case where we track every frequency line one by one, the covariance matrix of the associated Kalman filter has a size bigger than 100 × 100.
One can argue that recent high-end computers can handle the case where the size of the state vector is arbitrarily large. However, there may be a case where the ownship needs to track multiple targets. As we consume more power on tracking a single target, the total number of targets tracked by the ownship must decrease. Therefore, it is desirable to minimize the computational load for tracking a single target.
Moreover, 3D AFTMA has a large uncertainty in target range initially. For reducing the initial range estimate error in target-tracking problems, the RPEKF was applied [11,12]. The RPEKF uses multiple independent EKFs, each with a different initial range estimate, and the merged estimate is calculated as a weighted sum of individual EKF outputs. For addressing 3D AFTMA, we can apply multiple independent EKFs that are initialized at distinct range intervals, inspired by the RPEKF. However, as we consume more power on a single EKF track, we need to consume much more power as we use multiple EKF tracks. Therefore, it is desirable to minimize the computational load for a single EKF track.
Instead of tracking every frequency line one by one, we propose to estimate the average emitting frequency by setting the average frequency as the state vector in the filter. As the frequency measurements are averaged, the measurement noise decreases. Moreover, using the average frequency as our Kalman filter state, we can decrease the matrix size of the covariance matrix in the Kalman filter. This decreases the numerical error in matrix calculations, such as the calculation of an inverse matrix.
We show the outperformance of applying the average frequency using MATLAB simulations. As we use the average frequency line as our filter state, both the computational load and the RMSE decrease, compared to the case where we track every frequency line one by one.
At each sample-step k, one averages all frequency line measurements to obtain Then, (3) and (10) lead to where n f ,avg has a Gaussian distribution with mean 0 and standard deviation σ f ,avg = The standard deviation of n f ,avg is smaller than σ f . This implies that one can reduce the effect of noise as M k (the number of frequency lines) increases. In other words, as the frequency measurements are averaged, the measurement noise decreases.
Considering the 3D AFTMA problem with average frequency measurements in (15), the measurement equation at sample-step k is In (17), e k is the measurement noise with Gaussian distribution given as e k ∼ N(0, R k ).
Using (13) and (14), h(X k ) in (17) is calculated as Recall that in X k , X k [7] = f e,avg k , which appears in (3). In (18), f a (X k ) is the frequency measurement equation associated with (16) and is defined as whereṙ k is defined in (11).

Proposed Tracking Filter
We introduce the 3D AFTMA tracking filter applied in our paper. Recall that our tracking filter uses the relative state vector X k in (8). LetX k|k denote the estimation of X k based on all measurements up to sample-step k. In addition, let P k|k denote the error covariance ofX k|k . LetX k|k−1 denote the prediction of X k , based on all measurements up to sample-step k − 1. Let P k|k−1 denote the error covariance ofX k|k−1 .
Under the EKF with process model (9) and measurement model (17), one updates the state vectorX k|k and its covariance P k|k at every sample-step k. The detailed process of the EKF is discussed in [13].

A Frequency Line Appears or Disappears at Sample-Step k
Instead of estimating all emitting frequencies in S e k , this study runs the EKF for estimating f e,avg k , the average of emitting frequency. One presents a gating approach to determine whether a frequency average measurement f avg k in (16) can be applied for updating the EKF at sample-step k.
Recall that the sorted emitting frequency set S e k is defined in (2). Every element in S e k is sorted using (1). In practice, there may be a case where S e k−1 = S e k . This implies that a frequency line can appear or disappear at sample-step k.
Then, two gating conditions are as follows. and Here, max(S) denotes the maximum value of all elements in a set, say S. Satisfying both (22) and (23) indicates that we satisfy S e k−1 == S e k . In (23), β is a tuning parameter. In the case where (22) or (23) is not met at sample-step k, S e k−1 = S e k .

Reset the EKF Estimate for Emitting Frequency
If both (22) and (23) are satisfied at sample-step k, then one runs the EKF measurement update associated with f avg k . Otherwise, we reset the EKF estimate for emitting frequency, since frequency lines change at sample-step k.
In the case where the gating conditions are not satisfied at sample-step k, one can consider skipping the measurement update associated with f avg k . Skipping the measurement update is useful if the new frequency average measurement is associated with the existing frequency state X k [7]. However, our problem is that a new frequency line can suddenly appear at sample-step k and that the new frequency line may not be associated with the existing frequency state X k [7]. In the worst case, all new frequency lines appear at one sample-step, while the existing frequency lines disappear at the moment. This worst case scenario is simulated in Section 4.1.
Thus, we cannot apply skipping a measurement update in the case where the gating conditions are not satisfied at sample-step k. Instead, we reset the EKF estimate for emitting frequency, since frequency lines can suddenly appear at sample-step k.
One addresses how to reset the EKF estimate for emitting frequency, in the case where the gating conditions are not satisfied at sample-step k. Using (11) Here,ṙ(X k|k−1 ) indicates (11) when we useX k|k−1 instead of X k . Furthermore, (24) implies that this study uses the estimateX k|k−1 for the estimation of range rateṙ. It is acknowledged that the noise term n f ,avg in (16) is ignored in the derivation of (24). Recall that when using (16), n f ,avg has a Gaussian distribution with mean 0 and standard deviation σ f ,avg = Thus, σ f ,avg decreases as M k increases. This implies that as we use more frequency lines, the accuracy of (24) improves.
Recall that X k = (x k , y k , z k ,ẋ k ,ẏ k ,ż k , f e,avg k ) T denotes the state vector in our paper. The estimate vectorX k|k [1 : 6] is associated with the position and velocity of the target. One resetsX k|k [1 : 6] However, the error covariance for X k [7] needs to be reset. Using (24), one resets Cov(X k [7]) as (28) Here, σ f ,avg = σ f √ M k using (16). In addition, V max is the target's maximum speed that is assumed to be accessible in advance.

MATLAB Simulations
MATLAB simulations are presented to verify the effectiveness of the proposed 3D AFTMA filter. A simulation scenario lasts for 200 s. The sampling duration is set as dt = 0.1 s. According to [1], the sound speed in underwater environments is set as C = 1500 m/s. Let σ a define the standard deviation of a. In (7), one uses σ a x = σ a y = σ a z = 0.1 (m/s 2 ). In (17), the angle noise (σ b,1 or σ b,2 ) is set as 1 degree. In addition, the frequency noise σ f is set as 0.5 Hz. These noise parameters are feasible according to [1,2]. In a gating condition (23), we apply β = 100.
Let d 0 denote the relative distance between the ownship and the target at sample-step 0. Let d e = 10,000 (m) denote the range estimate error at sample-step 0. In the EKF,X 0|0 is initialized asX Here, b 1,0 and b 2,0 denote the azimuth and the elevation at sample-step 0, respectively. In addition, f avg 0 denotes the average frequency measurement at sample-step 0. See (16) for the equation of f avg k . The error covariance P 0|0 is set as Here, ).
The computational load is an evaluation indicator of algorithms. Let CompLoad denote the computational time for all MC simulations under MATLAB simulations.

Scenario 1
Scenario 1 is as follows. Initially, the ownship is positioned at (0, 0, 100) in meters. The ownship moves with a constant velocity (5, 0, 0) in m/s for 1000 sample-steps. From sample-step 1000 to 1150, the ownship performs the coordinates turn (CT) in the xy-plane with turn rate ω = 2 degrees per second. Here, the ownship's CT is modeled using Here, F CT is Here, S w = s(ω×dt) ω×dt , and C w = 1−c(ω×dt) ω×dt . Initially, the target is positioned at (2000, 0, 0) in meters. We consider a constantvelocity target, whose speed is set as Sp = 2 m/s. In (4), the target's velocity vector is set as Here, θ = π, and ψ = 3π 4 is used. Figure 1 depicts Scenario 1. In this figure, the target's trajectory is shown as red circles. In addition, the ownship's trajectory is shown as blue circles. A black asterisk denotes the initial position of the target or the ownship.  For comparison, one simulates the case where one does not reset the EKF estimate. In other words, one simulates the case where the filter reset strategy in Section 3.2 is not applied. Associated with Scenario 1 in Figure 1, NoReset in Figure 3 depicts RMSE k (33) in the case where the filter reset strategy in Section 3.2 is not applied. The RMSE suddenly increases at 1000 s, since frequency lines change at this moment. For NoRest, CompLoad is 33 s.
For another comparison, one simulates the ideal case where all frequency lines exist at all sample-steps. Associated with Scenario 1 in Figure 1,  Associated with the scenario in Figure 4, NoReset in Figure 5 depicts RMSE k in the case where the filter reset strategy in Section 3.2 is not applied. In NoReset, we use the average frequency line f avg k as our filter state. Considering NoReset, CompLoad is 34 s. Using (12) as the measurement equation, we can build a 3D AFTMA filter for tracking every frequency line one by one. For tracking every frequency line one by one, we initialize the EKF as follows. Since there are M frequency lines in total,X M 0|0 is initialized aŝ Associated with the scenario in Figure 4, useLarge in Figure 5 depicts RMSE k as we track every frequency line one by one. Considering useLarge, CompLoad is 64 s. We track every frequency line one by one, and the computational load increases compared to the case where we use the average frequency line f avg k as our filter state. Moreover, the RMSE increases compared to the case where we use the average frequency line as our filter state. Note that using the average frequency as our filter state, we can decrease the matrix size of the covariance matrix in the Kalman filter. This decreases the numerical error in the matrix calculations, such as the calculation of an inverse matrix. Figure 6 depicts Scenario 2. In this scenario, the trajectory of the ownship is identical to that in Scenario 1. A constant-velocity target starts from (500,0,0), and it moves with speed Sp = 5 m/s. In (4), the target's velocity vector is set as

Scenario 2
Here, θ = π 6 , and ψ = 3π 4 is used. In Figure 6, the target's trajectory is shown as red circles. In addition, the ownship's trajectory is shown as blue circles. A black asterisk denotes the initial position of the target or the ownship.
Associated with Scenario 2 in Figure 6, multiple frequency lines f m k (m ∈ {1, 2, . . . , M k }) are depicted with red color in Figure 7. The average frequency lines f avg k are depicted in blue. The frequency lines disappear and appear occasionally.  For comparison, one simulates the case where the filter reset strategy in Section 3.2 is not applied. Associated with Scenario 2 in Figure 6, NoReset in Figure 8 depicts RMSE k in the case where the filter reset strategy in Section 3.2 is not applied. The RMSE of NoReset suddenly increases at 500 s, since frequency lines change at this moment (see Figure 7). Considering NoReset, CompLoad is 34 s.
For another comparison, one simulates the ideal case where all frequency lines exist at all sample-steps. Associated with Scenario 2 in Figure 6, multiple frequency lines f m k (m ∈ {1, 2, . . . , M k }) are depicted with red color in Figure 9. The average frequency line f avg k is depicted in blue. No frequency lines disappear in this scenario.  Figure 10. RMSE of the scenario in Figure 9. NoReset depicts RMSE k in the case where the filter reset strategy in Section 3.2 is not applied. In NoReset, we use the average frequency line f avg k as our filter state. In this figure, useLarge depicts RMSE k as we track every frequency line one by one.
Using (12) as the measurement equation, we can build a 3D AFTMA filter for tracking every frequency line one by one. For tracking every frequency line one by one, we initialize the EKF state vectorX M 0|0 using (37). Associated with the scenario in Figure 9, useLarge in Figure 10 depicts RMSE k as we track every frequency line one by one. CompLoad under useLarge is 65 s. Since we track every frequency line one by one, the computational load increases compared to the case where we use the average frequency line as our filter state. Moreover, the RMSE of useLarge increases compared to the case where we use the average frequency line as our filter state.

Conclusions
In the 3D AFTMA problem, the ownship tracks a 3D target based on frequency, elevation, and azimuth measurements. This study addresses 3D AFTMA problems so that the ownship can track a target while measuring the target's sound with multiple frequency lines. Instead of tracking every frequency line one by one, we propose to estimate the average emitting frequency by setting the average frequency as the state vector in the filter. Moreover, our tracking filter handles the case where some frequency lines may disappear and appear occasionally. MATLAB simulations demonstrate the performance of the proposed tracking approaches. In the future, we will perform experiments utilizing a real ownship platform to demonstrate the proposed tracking approaches.