Tracking of Multiple Moving Sources Using Recursive EM Algorithm

We deal with recursive direction-of-arrival (DOA) estimation of multiple moving sources. Based on the recursive EM algorithm, we develop two recursive procedures to estimate the time-varying DOA parameter for narrowband signals. The ﬁrst procedure requires no prior knowledge about the source movement. The second procedure assumes that the motion of moving sources is described by a linear polynomial model. The proposed recursion updates the polynomial coe ﬃ cients when a new data arrives. The suggested approaches have two major advantages: simple implementation and easy extension to wideband signals. Numerical experiments show that both procedures provide excellent results in a slowly changing environment. When the DOA parameter changes fast or two source directions cross with each other, the procedure designed for a linear polynomial model has a better performance than the general procedure. Compared to the beamforming technique based on the same parameterization, our approach is computationally favorable and has a wider range of applications.


INTRODUCTION
The problem of estimating the direction of arrival (DOA) of plane waves impinging on a sensor array is of fundamental importance in many applications such as radar, sonar, geophysics, and wireless communication.The maximum likelihood (ML) method is known to have excellent statistical performance and is robust against coherent signals and small sample numbers [1].However, the high computational cost associated with ML method makes it less attractive in practice.
To improve the computational efficiency of the ML approach, numerical methods such as the expectation and maximization (EM) algorithm [2] were suggested in [3,4,5].Recursive procedures based on the recursive EM (REM) algorithm for estimating constant DOA parameters were discussed in [6,7].Similar procedures for tracking multiple moving sources were studied in [8,9].In [9], the authors focused on narrowband sources and assumed known signal waveforms.
The REM algorithm is a stochastic approximation procedure for finding ML estimates (MLEs).It was first suggested by Titterington [10] and extended to the multidimensional case in [6].As it was pointed out by Titterington, REM can be seen as a sequential approximation of the EM algorithm.The gain matrix of REM is the inversion of the augmented data information matrix.Through proper design of the augmentation scheme, the augmented data and the corresponding information matrixes usually have a simplestructure [2].In this case, the REM algorithm is very easy to implement.For constant parameters, estimates generated by REM are strongly consistent and asymptotically normally distributed.For time-varying parameters, the tracking ability of a stochastic approximation procedure depends mainly on the dynamics of the true parameter, gain matrix, and step size [11].
Based on REM, we will derive two recursive procedures for estimating time-varying DOA.The first procedure does not require any prior knowledge on the motion model.The only assumption is that the unknown parameter changes slowly with time.The second procedure assumes that the time-varying DOA parameter θ(t) is described by a linear polynomial of time.This model is important since a smooth function θ(t) can be approximated by a local linear polynomial in a short-time interval [12].The procedure reported in [8] employs a decreasing step size to estimate the polynomial coefficients.However, since the DOA parameter θ(t) and the log-likelihood function change with time, a decreasing step size may not capture the nonstationary feature of the underlying system over a long period.To overcome this problem, we suggest a constant step size to be used in the algorithm.It is noteworthy that both procedures are aimed at maximizing the expected concentrated likelihood function [13].Introducing a linear polynomial model implies increasing the dimension of the parameter space.With the additional degree of freedom, the procedure designed for a linear polynomial model should perform better than the general one.
In contrast to methods based on subspace tracking [14] or two-dimensional beamforming [12], our approach can be easily generalized to wideband cases including underwater acoustic signals.Unlike the Kalman-type algorithms [15], the recursive procedures considered here have a much simpler implementation.
This paper is outlined as follows.We describe the signal model and the REM algorithm briefly in Sections 2 and 3. Section 4 presents two recursive procedures for localizing moving sources.Simulation results are discussed in Section 5. We give concluding remarks in Section 6.

PROBLEM FORMULATION
Consider an array of N sensors receiving M far-field waves from unknown time-varying directions where the steering matrix To avoid ambiguity, we assume that M < N. The signal waveform T ∈ C M×1 is considered unknown and deterministic.(•) T denotes the transpose of a vector.Furthermore, the noise process u(t) ∈ C N×1 is independent identically complex and normally distributed with zero mean and covariance matrix νI, where ν represents the unknown noise spectral parameter and I is the identity matrix.
In the following, we assume that the number of sources M is known.Standard procedures based on minimum description length (MDL) criteria [16] or multiple hypothesis testing [7] can be used to determine M. The problem of interest is to estimate the time-varying DOA parameter θ(t) recursively from the observation x(t).We assume that a good initial estimate θ 0 is available at the beginning of the recursion.

RECURSIVE PARAMETER ESTIMATION USING INCOMPLETE DATA
The REM algorithm suggested by Titterington is a stochastic approximation procedure for finding MLEs.As pointed out in [10], there is a strong relationship between this procedure and the EM algorithm [2].Using Taylor expansion, Titterington showed that, approximately, REM maximizes EM's augmented log likelihood sequentially.The unknown parameter is considered as constant in [10].In the fixed parameter case, a properly chosen decreasing step size is necessary to ensure strong consistency and asymptotic normality of the algorithm [10,17].Suppose x(1), x(2), . . .are independent observations, each with underlying probability density function (pdf) f (x; ϑ), where ϑ denotes an unknown constant parameter.The augmented data associated with the EM algorithm y(1), y(2), . . . is characterized by the pdf f (y; ϑ).According to [2], the augmented data y(t) is so specified that M(y(t)) = x(t) is a many-to-one mapping.Let ϑ t denote the estimate after t observations.The following procedure is aimed at finding the true parameter ϑ which may coincide with the MLE in the asymptotic sense [18]: where t is a decreasing step size and represent the augmented information matrix and gradient vector, respectively.∇ ϑ is a column gradient operator with respect to ϑ.We assume that both (4) and (5) exist.Under mild conditions, the estimates generated by (3) are strongly consistent, asymptotically normally distributed.In view of the well-known singularities and multiple maxima that are on likelihood surfaces, one could not of course expect consistency irrespective of the starting point [10].
The augmented data y usually has a simpler structure than the observed data x.Therefore, the augmented data information matrix I EM (ϑ t ) is easier to compute and invert than the observed data information matrix Although REM does not have the optimal convergence rate in the Cramér-Rao sense as the following procedure [10]: it is much easier to implement than (6).Using I EM (ϑ t ) −1 as the gain matrix is a tradeoff between the convergence rate and computational cost.
When the parameter of interest varies with time, a decreasing step size such as t = t −α , 1/2 < α ≤ 1, cannot capture the nonstationary feature of the underlying system.A classical way to overcome this difficulty is to replace t with a constant step size .In general, a large step size reduces the bias and increases the variance of the estimates [11].
A small step size has opposite effects.Since the time-varying parameter ϑ(t) may follow a complicated dynamics, an exact investigation of the convergence behavior of the algorithm is only possible when certain assumptions are made on the parameter model.More discussion about convergence properties of a stochastic approximation procedure in a nonstationary environment can be found in [11].

LOCALIZATION OF MOVING SOURCES
The REM algorithm with constant step size ( 7) is applied to estimate the time-varying DOA parameter θ(t).We start with a general case in which θ(t) changes slowly with time and then consider a linear polynomial model.

General case
From the signal model in Section 2, we know that the array observation x(t) is complex and normally distributed with the log likelihood function where ϑ = [θ(t) T s(t) T ν] T and (•) H denotes the Hermitian transpose.
According to (7), all elements in ϑ should be updated simultaneously.Since we are mainly interested in the DOA parameter θ(t) and since including {s(t), ν} in the recursion will complicate the gain matrix I EM (ϑ t ) −1 , the procedure (7) is only applied to θ(t).The estimate for signal waveform and noise level, denoted by T and ν t , respectively, is updated by computing their MLEs once the current DOA estimate is available.For simplicity, we use θ instead of θ(t) in the following discussion.
Taking the first derivative on the right-hand side of ( 8) with respect to θ m , we obtain the mth element of the gradient vector γ(x(t), ϑ t ) [17]: where The augmented data y(t) is obtained by decomposing the array output into its signal and noise parts.Formally it is expressed as The augmented data associated with the mth signal (1) Calculate the gradient vector γ(x(t), θ t ) by ( 9) and the matrix I EM (θ t ) by ( 14). ( 2) Update DOA parameters by (3) Update the signal and noise parameters s t , ν t by (15).
is complex and normally distributed with mean d(θ m )s m (t) and covariance matrix Since the signals are decoupled through the augmentation scheme (10), I EM (ϑ t ) is an M × M diagonal matrix when we only consider the DOA parameter θ.By definition (4), the mth diagonal element of I EM (ϑ t ) is the conditional expectation of the second derivative of the augmented log likelihood which is given by where Once the estimate θ t+1 is available, the signal and noise parameters are obtained by computing their MLEs at current θ t+1 and x(t) as follows: where H(θ t+1 ) # is the generalized left inverse of the matrix H(θ t+1 ), P(θ t+1 ) ⊥ = I − P(θ t+1 ) is the orthogonal complement of the projection matrix P(θ t+1 ) = H(θ t+1 )H(θ t+1 ) # , and C x (t) = x(t)x(t) H .Given a constant step size , the number of sources M, and the current estimate θ t , the (t + 1)st recursion of the algorithm proceeds as shown in Algorithm 1.
The REM algorithm is applied to estimate θ 0 and θ 1 .For notational simplicity, we define the extended DOA parameter as Similarly to the procedure presented in Section 4.1, REM is only applied to update the DOA parameter Θ rather than Based on this approach, the 2mth and (2m + 1)st element of the gradient vector γ(x(t), ϑ t ) are given by respectively, where d (Θ t m ) = ∂d(θ m )/∂θ m | θm=θ t 0m +tθ t 1m .Note that θ is calculated at the current estimate Θ t according to the linear model (16).
Because each source is described by two unknown parameters, the augmented data information matrix becomes block diagonal.Unfortunately, this matrix is singular under current parameterization.To avoid singularity and simplify the recursion, rather than using this block diagonal matrix in the recursion directly, we consider an alternative matrix ĨEM (ϑ t ) which is the diagonal part of I EM (ϑ t ). Let According to the augmentation scheme specified above, the 2mth and (2m + 1)st diagonal components of ĨEM (Θ t ) are given by respectively.
Similarly to the general case, the signal and noise parameters are updated by (15) once the estimate Θ t+1 is available.The parameter θ t+1 in ( 15) is replaced by Θ t+1 .
Given the step size , the number of sources M, and the current estimate Θ t , the (t + 1)st recursion of the algorithm proceeds as shown in Algorithm 2.
For simplicity, the REM for the general case and the REM for the linear polynomial model are referred to as "REM I" and "REM II," respectively.
From ( 9), (14), and (15), the computational complexity of REM I lies approximately between O(MN + MN 2 ) and O(MN + N 3 ).The dominant term MN 2 (or N 3 ) is associated with s t+1 given by ( 15) which is a solution to a least square (LS) problem.Different LS algorithms yield different computational loads [19].Due to the increased number of unknowns, REM II requires twice as many computations as REM I in computing the gradient vector and augmented information matrix.Clearly, REM II is computationally more efficient than the localpolynomial-approximation (LPA) based beamforming technique [12] whose computational complexity is given by O(NTLP) where T represents the number of snapshots, L denotes the number of points in the angular search domain, and P denotes the number of angular velocity search domain.
It was pointed out in [13] that REM for constant DOA estimation is indeed a recursive procedure for finding the maximum of the expected concentrated likelihood function where The constant step size considered in REM I captures the time-varying character of the likelihood function.Similarly, REM II is aimed at finding the maximum of L(θ).Using a different parameterization, such as a linear polynomial model implies increasing the dimension of the parameter space.With the additional degree of freedom, REM II is expected to have a better tracking ability.Later in Section 5 we will show that in critical situations where two source directions cross with each other, REM II provides more accurate estimates than REM I. Choosing a proper step size plays an important role in the algorithms' tracking ability.The optimal step size depends on the dynamics of the true parameters, for instance, rate of change.Interested readers can find general guidelines in [11] and an adaptive procedure designed for REM with a decreasing step size in [20].

Extension to broadband signals
The algorithms presented previously are derived under the narrowband signal assumption.Extension to the broadband case is straightforward.From the asymptotic theory of Fourier transform [21], we know that each frequency bin is asymptotically independent of each other [22].The log likelihood function associated with the broadband signal is a sum of the log likelihoods of individual frequency bins.Correspondingly, the gradient vector and augmented information matrix can be easily obtained by adding up the gradient vectors and augmented data information matrices of relevant frequency bins.Similarly to the narrowband case, the signal and noise parameters at each frequency are updated by calculating their MLEs once the current DOA estimate is available.

SIMULATION
The proposed algorithms are tested by numerical experiments.In the first part, we consider REM algorithms' application in narrowband and broadband cases.In the second part, we compare REM II with the LPA-based beamforming technique [12].

Recursive EM algorithms I and II
The narrowband signals generated by three sources of equal power are received by a uniformly linear array of 15 sensors with interelement spacings of half a wavelength.The signalto-noise ratio (SNR), defined as 10 log(s m (t) 2 /ν), m = 1, 2, 3, is kept at 10, 20 dB.The motion of the moving sources is described by the linear polynomial model (16).Three different parameter sets {θ 0 , θ 1 } are assumed in the experiments.Each experiment performs 200 trials.
In the first experiment, we consider relatively fast moving sources.The true parameters are given by θ where θ 1 is mea- sured by degrees per time unit.In order to get a good insight into the tracking behavior, the same initial values are used in all trials.We applied LPA-based beamforming to 20 snapshots to obtain the initial estimates θ 0 0 = [10.5 • , 59.5 The initial estimate for REM I is given by θ 0 0 .Both algorithms use a constant step size = 0.6.Figures 1 and 2 present the true values of θ and an example of estimated trajectories.As shown in both figures, two source directions cross with each other at t = 32.Obviously, the recursive procedure designed for the most general case cannot follow fast moving sources at all.In contrast, the estimated trajectory obtained by REM II is very close to the true one.Figures 3 and 4 show the root mean square errors (RMSEs) of the DOA estimates, defined as 2 , averaged over 200 trials.Since REM I fails to track the moving sources, the corresponding RMSE grows with increasing time.On the other hand, the RMSE associated with REM II decreases slightly at the beginning of the recursion and then remains almost constant.Comparing Figures 3 and 4, one can observe that SNR = 20 dB has a slightly lower RMSE than SNR = 10 dB.
The second experiment involves three slowly moving sources.The true parameter values are given by θ Note that the angular velocity θ 1 is approximately 1/10 of that considered in the previous experiment.We applied the ML method to obtain the initial estimates θ 0 0 = [30.1 • , 50.8 • , 60.9 • ].Because the angular velocity is very small compared to that in the previous experiment, we take θ 0 1 = [0 • , 0 • , 0 • ] as the initial value for θ 1 .The initial estimate for REM I is given by θ 0 0 .Both algorithms use a constant step size = 0.6.Figures 5  and 6 present the true and estimated trajectories obtained by REM I and REM II.Similarly to the first experiment, two source directions cross with each other at t = 126.The estimated trajectory by REM I is close to the true one when no  crossing happens.Between t = 100 and t = 230, where two source directions cross with each other, the estimated trajectories associated with the first two sources do not get close to each other.Instead, they just depart in the vicinity of t = 126.For the same scenario, REM II provides a more accurate estimate.Figure 6 shows that the crossing point causes a larger deviation from the true trajectory.Due to a higher sensitivity to the variation of angular velocity at the crossing point, the estimated trajectory in Figure 6 is slightly worse than that in Figure 2. Comparison of Figures 7 and 8 with Figures 3 and  4 shows an overall lower RMSE in this scenario.Although REM I provides more reliable estimates than in the first experiment, REM II still outperforms REM I.In the third experiment, three sources move slowly with different speeds but do not cross with each other.The true parameters are given by θ We use a constant step size = 0.6.Both algorithms have good tracking ability.Figures 9 and 10 show that RMSE is the lowest among all three scenarios.REM II has a better performance than REM I.While REM II has a better performance at higher SNR, REM I seems to be less sensitive to SNRs in all three scenarios.
In addition to the narrowband signals, we also applied REM I and REM II to broadband signals with 3 frequency  bins.The scenario similar to the second experiment leads to the results presented in Figures 11 and 12.The estimates behave similarly to the narrowband case.Comparison of RM-SEs shows that using more frequency bins leads to higher accuracy.

Comparison with LPA beamforming
We compare REM II with the LPA-based beamforming approach suggested by Katkovnik and Gershman [12].Both algorithms assume the motion model (16).In the first experiment, the narrowband signals are generated by the following parameter sets    RMSE associated with REM II changes slowly over time.While estimates of θ 0 remain constant, the estimates of θ 1 become more accurate with increasing recursions.Also, we can observe that while LPA beamforming provides an overall better θ 0 estimates and better angular velocity estimates at beginning of the recursion, REM II improves θ 1 estimates with increasing time and has less fluctuations.
Compared with the Cramér-Rao bounds (CRBs) [23], one realizes that an REM II is certainly not an efficient estimator.However, the ML approach suggested in [23], whose estimation accuracy is close to CRB, is a batch processing and requires a complicated multidimensional search procedure.In the second experiment, REM II provides much more accurate estimates than LPA beamforming.Figure 17 shows that LPA beamforming even fails to follow the moving sources.We can observe in Figure 18 that REM II has lower RMSE in both θ 0 and θ 1 estimation.Consequently, as shown in Figures 19 and 20 the resulting DOA estimates are much better than LPA beamforming.In both experiments, the computational time needed for LPA beamforming is about 800 times as high as that required by REM II due to the twodimensional search procedure.
We conclude that REM I is suitable for tracking slowly time-varying DOA parameters, REM II performs well for both slowly and fast moving sources.Both procedures generate accurate estimates when there is no crossing point.When two source directions coincide with each other, the steering matrix H(θ) becomes rank deficient.The signal waveform s(t) cannot be determined properly.Consequently the DOA parameter cannot be estimated accurately.In this case, regularization is needed [23].Since REM II incorporates a linear polynomial model, it has a better tracking ability than REM I when this critical situation occurs.Compared to LPA beamforming, our method has a clear computational advantage.It provides comparable results with LPA beamforming in the fast moving sources case and outperforms LPA beamforming in the slow moving source case.In addition, REM is applicable to both narrowband and broadband signals.

CONCLUSION
We addressed the problem of tracking multiple moving sources.Two recursive procedures are proposed to estimate the time-varying DOA parameter.We applied the recursive EM algorithm to a general case in which the motion of the sources is arbitrary and a specific case in which the motion of sources is described by a linear polynomial model.Because of the simple structure of the gain matrix, the suggested procedures are easy to implement.Furthermore, extension of our approaches to broadband signals is straightforward.Numerical experiments showed that our approaches provide excellent results in a slowly changing environment.When the DOA parameter changes fast or two source directions cross with each other, the procedure derived for a linear

Figure 2 :
Figure 2: True trajectory (solid lines) and estimated trajectory (nonsolid lines) by REM II for the narrowband case.SNR = 20 dB.

Figure 12 :
Figure 12: RMSE of θ versus time for the broadband case.SNR = 10 dB.