A Fast ML-Based Single-Step Localization Method Using EM Algorithm Based on Time Delay and Doppler Shift for a Far-Field Scenario

This study discusses the localization problem based on time delay and Doppler shift for a far-field scenario. The conventional location methods employ two steps that first extract intermediate parameters from the received signals and then determine the source position from the measured parameters. As opposed to the traditional two-step methods, the direct position determination (DPD) methods accomplish the localization in a single step without computing intermediate parameters. However, the DPD cost function often remains non-convex, thereby it will cost a high amount of computational resources to find the estimated position through traversal search. Weiss proposed a DPD estimator to mitigate the computational complexity via eigenvalue decomposition. Unfortunately, when the computational resources are rather limited, Weiss’s method fails to satisfy the timeliness. To solve this problem, this paper develops a DPD estimator using expectation maximization (EM) algorithm based on time delay and Doppler shift. The proposed method starts from choosing the transmitter-receiver range vector as the hidden variable. Then, the cost function is separated and simplified via the hidden variable, accomplishing the transformation from the high dimensional nonlinear search problem into a few one dimensional search subproblems. Finally, the expressions of EM repetition are obtained through Laplace approximation. In addition, we derive the Cramér–Rao bound to evaluate the best localization performance in this paper. Simulation results confirm that, on the basis of guaranteeing high accuracy, the proposed algorithm makes a good compromise in localization performance and computational complexity.


Introduction
The transmitter localization is a classic problem in wireless communication systems. It is well known that the conventional location approaches are composed of two separate steps: (1) The intermediate parameters are estimated through the signals. (2) The source position is determined from the measured parameters. In the past three decades, a great deal of work has been done in this field. Generally, the intermediate parameters (e.g., angle of arrival (AOA) [1], time of arrival (TOA) [2], time difference of arrival (TDOA) [3], Doppler shift [4][5][6] and frequency difference of arrival (FDOA) [7]) are usually estimated by the maximum likelihood (ML)-based method [8] or the subspace data fusion (SDF)-based method [5,9], and the position is mainly determined by an iterative algorithm [10] or a closed-form algorithm [11]. It is worth mentioning that Doppler information is a key parameter in the location problem based on the moving receiver. Gajewski [5] uses SDF criterion based on the Doppler effect to locate multiple emission sources. Since multipath propagation and Doppler effect are dominant factors to deteriorate localization performance, Kelner [6] uses the signal Doppler frequency method to resist non-line-of-sight (NLOS) conditions. As for the conventional approaches, it must be emphasized that the intermediate parameters are extracted by ignoring the constraint that all measurements must correspond to the same geolocation of the emitter, and more errors will be introduced in two separate processes of the conventional methods [12]. As a result, the conventional two-step methods cannot yield satisfactory location accuracy. For this reason, direct position determination (DPD) techniques that exploit the intrinsic constraint and determine the source position from the received signals directly were developed.
The DPD algorithms have been intensively investigated in recent years. An important classification of localization methods is related to propagation conditions, and there are two general cases: (1) multipath propagations with LOS and NLOS; (2) an additive white Gaussian noise (AWGN). In the first case, the methods [13][14][15][16] for LOS or NLOS environments have been deeply developed. Note that Du [16] proposed a DPD estimator for a novel localization architecture in multipath propagations, called the "Multiple Transponders and Multiple Receivers for Multiple Emitters Positioning System". The second case is the most common condition, and the following discussion is also conducted under this assumption. Weiss [17] first proposed an SDF-based DPD algorithm with the utilization of orthogonality between signal subspace and noise subspace to estimate the emitter position. Chen [18] developed a multi-target DPD approach using subspace based on compressive sensing, reaching a high probability of locating the emitter without knowing the number of targets. Even if low computational complexity appears in the SDF-based DPD method, localization performance is deteriorated in low signal-to-noise ratios (SNRs), failing to reach the corresponding Cramér-Rao bound (CRB). To enhance the performance, ML-based DPD approaches were developed. Since the DPD cost function is often non-convex, traversal search is required to find the extremum. However, nuisance parameters (e.g., unknown transmission time [19] and timing errors [20]) will result in high dimensional search, which is impractical in a real-time application. To mitigate the computational load of traversal search, Weiss [21] proposed an ML-based approach via finding the maximum eigenvalue of the Hermitian matrix associated with the cost function, which also exhibits high localization accuracy. It must be emphasized that Weiss's method is an ideal solution to account for localization accuracy and computation complexity in existing DPD algorithms. Moreover, the location performance can be further enhanced by utilizing signal properties. These DPD algorithms [22][23][24] can obtain superior localization precision by exploiting the constant modulus, orthogonal frequency division multiplexing, and the cyclostationary properties of signals. Unfortunately, the above DPD methods are not fast enough in the presence of limited computing resources, which is often a reality in moving receivers (i.e., airplanes or unmanned aerial vehicles (UAVs)). This limitation results in non-timeliness performance for moving emitter localization. On the other hand, DPD methods based on the ML criterion guarantee superior localization accuracy. Therefore, an ML-based DPD method with rapidity needs to be studied.
The expectation maximization (EM) algorithm is an attractive method of estimating the ML result when data can be divided into "incomplete data" and "complete data" in the model. In the past three decades, the EM algorithm has provided an excellent way to solve machine learning problems (i.e., speech processing and recognition [25] and image identification and restoration [26]). Via choosing the appropriate hidden variable, the EM algorithm decouples a high-dimensional search problem into a few subproblems with much lower computational complexity. In sight of the advantages of the EM algorithm, many scholars have applied it to the parameter estimation domain. Mada [27] uses the EM algorithm to solve the multi-source localization problem, leading to a soft computational load. Moreover, Lu estimates the source position via the EM algorithm for spatially nonwhite noise [28] and nonuniform noise [29], respectively, further demonstrating the effectiveness of the EM model in harsh scenarios. However, the DPD methods have not received the same treatment on the EM application, and only Tzoreff [30] discussed a DPD method based on the EM algorithm. Unfortunately, this method cannot suit for the moving-receiver scenario. As a result, there is a great demand for developing a DPD method using the EM algorithm for moving receivers.
This paper proposes a fast ML-based DPD algorithm using the EM algortihm based on time delay and Doppler shift, and the main processing steps are exhibited as follows: (1) The transmitter-receiver range vector is selected as the hidden variable, successfully leading the separation and simplification of the ML cost function. (2) With the help of Laplace approximation, the high-dimensional multi-parameter search of the prescribed ML estimator is decoupled into a closed form of the transmitter position and a line search of transmitter-receiver distance as well as transmitted time. Therefore, the expressions of EM repetition is determined. (3) Iteration of the EM expressions, which alternately updates parameters in E-step and M-step, is required until the norm of the difference between the adjacent estimated position converges to a user's predefined number.
In summary, the main contribution of this paper is the improvement of the prescribed ML DPD estimator via the EM algorithm in moving-receiver application for a far-field scenario, leading to a high effectiveness to find the global extreme. In addition, we also have derived the CRB to exhibit the best localization performance in theory. The rest of this paper is organized as follows. Section 2 lists the notations used in this paper. Section 3 constructs the signal model and formulates the problem. Section 4 discusses the DPD methods, including the previous method and the proposed method, and then makes a computational complexity comparison among different methods. Section 5 presents several numerical simulations to verify the theoretical analysis. Finally, Section 6 draws the conclusions.

Notations
In this section, some mathematical notation explanations that will be used through this paper are listed in Table 1.

Signal Model
This paper considers that L moving receivers intercepts the transmitted signal, and the signal is partitioned into K short intervals. Note that the antennas of the emitter and receivers are omnidirectional, and the positions of the emitter and receivers are determined in a two-dimensional coordinate system. In order to fully introduce the DPD model, three assumptions are included as follows: Assumption 1. The receivers are moving. Let p l,k and v l,k denote the position and velocity of the lth receivers at the kth interception interval, which are precisely known to us. They are constant at each interception interval. Assumption 2. The stationary emitter, denoted by p fi px, yq, locates in the far-field of the moving receivers. Assumption 3. The propagation channel is an AWGN channel, and time delay as well as Doppler information are used in the signal model.
Based on the above assumptions, the complex signals y l,k ptq observed by the lth receiver at the kth interception interval at time t is expressed as [31] y l,k ptq " b l,k s k`t´τ l,k ppq´t 0˘e where t 0 is the signal transmission time, and during the kth interception interval, s k ptq is the complex envelope of the emitter, b l,k denotes the channel attenuation, n l,k ptq is the white Gaussian, τ l,k ppq is the propagation time between the emitter and the lth receiver, and Doppler frequency f l,k observed between the emitter and the lth receiver is given by [21] f l,k fi f c`fc µ l,k ppq where f c is the nominal frequency of the transmitted signal, assumed known, with where c is the radio wave propagation speed. It is assume that each receiver performs a down conversion of the intercepted signal by the nominal frequency. Thus, after down conversion the Doppler frequency can be replaced by After sampling at t " nT s , the received signal can be shown as y l,k pnq " b l,ks k`n T s´τl,k ppq´t k˘e j2π f l,k nT s`ñ l,k pnq n " 1, 2, . . . , N for l " 1,¨¨¨, L, and k " 1,¨¨¨, K, where N denotes the number of sample points at each interval. Then, Equation (5) can be expressed by a vector form as withÑ " r1, 2, . . . , Ns T . We defines k " rs k pT s q ,s k p2T s q ,¨¨¨,s k pNT s qs T , and the Fourier coefficients ofs τ l,k ands k satisfys F is the discrete Fourier transform operator. Therefore, Equation (6) is rewritten as where

Direct Position Determination Methods
This section discuss the DPD methods, which locate the emitter directly through the raw data. We first discuss the ML-based DPD method requiring traversal search, which can receive the best location accuracy. Weiss [19] proposed a fast ML-based DPD method based on eigenvalue decomposition. Unfortunately, the above methods cannot work effectively in the presence of limited computational sources. Thus, in light of the EM idea in estimation theory [27], a fast DPD method using EM algorithm will be proposed.

Previous Work
The DPD optimization based on time delay and Doppler shift is established, which was first introduced by [21]. The received signalỹ l,k is a complex Gaussian random vector. Hence, the likelihood function forỹ can be formulated by [21] l pθq " 1 ( 11) where σ 2 denotes the noise power, and θ " " t 0 , b T , p T ‰ T denotes all unknown parameters, The associated logarithmic likelihood function can be written as Therefore, the estimation of noise power σ 2 iŝ By substituting Equation (11) into Equation (10), the estimation of parameter θ can be determined bŷ with f pθq " Since a high-dimensional nonlinear problem appears in Equation (14), it is difficult to compute the closed-form solution ofθ. Thus, traversal searches are required among these stray parameters to obtain the best performance. However, this technique will take a long time to find the extreme corresponding to the emitter position, which is not efficient in practical application.

EM Algorithm Review
The EM algorithm is an approach to iterative computation of the ML problem when the observations are regarded as incomplete data. In each iteration, it includes an expectation step and a maximization step. The meaning "incomplete data" reveals that there are two kinds of data: the incomplete data Y and the complete data X. More specifically, Y is the observed data, and X is the corresponding hidden data, (not observed directly). We denote the estimated parameter θ and an expression Y " f pXq, which shows a many-one mapping from X to Y. Via the Bayesian rule, we have [30] L pθq fi log p Y pY; θq " log p X,Y pX, Y; θq´log p X{Y pX{Y; θq (16) where L pθq is the logarithmic likelihood function. Since p Y pY; θq is independent to X, the conditional expectation operation with associated with p X{Y`X {Y; θ 1˘w ill not make change to L pθq. The conditional expectation of Equation (14) can be expressed by where E " E pX{Y;θ 1 q t‚u and θ 1 denotes an arbitrary value of θ.
We define Thus, Equation (17) can be rewritten as Based on the ML criterion, we can maximize L pθq to estimate the unknown parameters. Since p X{Y pX{Y; θq is generally unknown to us, we need to approximate W´θ, θ 1¯.
Both sides of Equation (19) subtract L´θ 1¯y ielding With the help of Gibbs inequality, we find W´θ, θ 1¯ě W´θ 1 , θ 1¯. Thus, we obtain The above expression shows that an increment of Q associated with θ also ensures an increment of L. Therefore, the maximization problem of L is transformed into the maximization problem of Q. The EM algorithm contains two steps, and each iteration process can be given by An iteration cycle exists in Equation (20). We obtain θ pi`1q in a current iteration step, and it will be the initial value to repeat the EM operation of Equation (20) in the next iteration step. When Q converges, the iteration stops.

Derivation of the EM-DPD Algorithm
After the above analysis, we choose the received signalsỸ " ı T as the hidden data, respectively. We assume thatx l,k is a Gaussian vector, and the probability distribution function (PDF) ofX can be shown as pX`X; p˘" ( 23) where ε l,k ppq " p´p l,k and σ 2 x is the variance of x. We define G l,k fi C l,ksk , and Equation (6) is rewritten asỹ We assume thatñ l,k are complex Gaussian vectors with mean 0 N and covariance σ 2 I N . Therefore, the pdf ofỹ l,k is given by Assuming the signals observed by different receivers and different observations intervals are independent, we have pỸ`Ỹ; p, t 0˘" Note that the source coordinates p is the direct relation parameter to the observationsỸ. We need to separate this variable to the subsequent derivation. By introducingX in Equation (23) and defining φ " (26) can be rewritten as pỸ`Ỹ{X; φ˘" where G l,k " G l,k`xl,k , t 0˘f or brevity. From the analysis in Section 4.2.1, we can estimate the emitter position via maximizing the auxiliary function Q´θ, θ 1¯. Next we introduce the key procedure of the derivation.

Proposition 1.
The auxiliary function Q´θ, θ 1¯i s separated into Proof of Proposition. . Via the Bayesian rule, the joint probability pX ,Ỹ`X ,Ỹ; θ˘is expressed by the product of Equations (23) and (26), which depend only on p and only on φ, respectively. Since the logp‚q of pX ,Ỹ`X ,Ỹ; θ˘exists in Equation (18), the separation of Q´θ, θ 1¯c an be shown as By maximizing Equation (29), we obtain whereḠ l,k " EtG l,k u and E s " › › G l,k › › 2 . Substituting Equation (31) into Equation (29) yields where Ψ l,k " I N´E´1 sḠl,kḠ H l,k . We continue to maximize the above expression, yieldinĝ Similarly, maximizing Equation (28) After eliminating b in the optimization process, we obtain the ML estimated expression of t 0 and the closed-form solution of p through the EM procedure. Note thatḠ l,k in Equation (31) andx l,k in Equation (32) are need to estimate. This work is promoted by Laplace approximation, which is shown in Appendix A. We define r l,k fi M-step: For easy understanding, the main steps of the proposed method are exhibited in Algorithm 1.

Algorithm 1:
The main steps of the proposed method.

Input:
The observed data:ỹ l,k , the position, and the velocity of receiver: p l,k and v l,k , pl " 1,¨¨¨, L, k " 1,¨¨¨, Kq; 1. Choose a small positive number ε ą 0, and set the iteration counter to i " 1; 2. Set i=0, initialize p piq , t

Computational Complexity Analysis
Based on the above analysis, we have achieved the transformation from high-dimensional multi-parameter nonlinear problem into several optimization subproblems, whose computational complexity is substantially reduced. Nonetheless, the computational complexity differs among the traversal search method, Weiss's method [21], and the proposed method. The traversal search method requires a three-dimensional search to find the extreme of the cost function corresponding to the emitter position. Weiss's method mainly has a two-dimensional search of the maximum eigenvalue of a Hermitian matrix. The proposed method is dominated by the line searches of rl ,k pl " 1,¨¨¨, L, k " 1,¨¨¨, Kq in Equation (35) and t 0 in Equation (37).
To better exhibit the computational complexity, Table 2 lists the computational complexity of the traversal search method, Weiss's method, and the proposed method. Note that N p is the number of grid search points with respect to a line search, and I iter is the number of iterations of the proposed method. It must be emphasized that the key contributor to computational complexity is N p . Compared with the two-or three-dimensional searches in other methods, the proposed method only has a one-dimensional search, which reduces the computational complexity significantly.

. Numerical Experiments
In this section, several numerical experiments are reported to corroborate the theoretical analysis. All simulations results are based on 200 Monte Carlo trials. In this scenario, the receivers equipped with only one sensor are included, and the source emits a Gaussian random signal with a center frequency of f c = 200 MHz. The channel attenuation coefficient amplitude obeys a normal distribution with a mean of 1 and a standard deviation of 0.1, and the channel phase obeys a uniform distribution over r´π, πs. Unless otherwise specified, we collect N " 32 sample points in each interval at a sampling rate of f s = 400 kHz, use L = 3 receivers, perform a total of K = 5 observations, and set the velocity of receiver v = 300 m/s.
To examine the performance comparisons, we take simulations works with the following four algorithms 1. the proposed method in this study; 2. the traversal search method; 3. Weiss's method; 4. the TOA/Doppler two-step algorithm.
The details of the comparison algorithms are shown in Table 3. Additionally, the CRB presented in Appendix B is also included in this part, providing a theoretical best performance benchmark. Moreover, root mean square error (RMSE) and cumulative distribution function (CDF) are adopted to evaluate localization accuracy in this paper.
We now examine the localization performance of the proposed algorithm for three different scenarios. The target locates at [5,4] km, and the receivers move along different trajectories, which can be found in Figure 1. As shown in Figure 2, the performance of our algorithm is not sensitive to the receiver trajectories because there is no significant difference between the RMSE curves of our algorithm for different scenarios. Additionally, as expected, the proposed method for Scenario (a) has better localization precision in low SNRs. On the other hand, although the RMSE of our algorithm is far from that of the corresponding CRB in low SNRs, it has a substantial advantage in computation time (see the descriptions for Table 4).  It must be emphasized that the next simulations are all based on Scenario (a) in Figure 1. Firstly, we continue to conduct the simulation for algorithm performance comparison versus SNRs. The RMSEs of all algorithms can be easily found in Figure 3. Unsurprisingly, all DPD methods outperform the two-step method, and the traversal search method shows the best localization performance especially in low SNRs. Although the performance of the proposed method is slightly inferior than that of Weiss's method in low SNRs, it is more efficient. However, the other two DPD methods will cost a high amount of time in the presence of limited computing resources, losing its timeliness to moving target. The detailed runtime information can be seen in Table 4. Secondly, to better exhibit the computational complexity of each algorithm, the runtime is investigated. From the results in Table 4, we observe that the two-step method requires the least runtime, but its localization performance is deteriorated in low SNRs (see Figure 3). Due to the high-dimensional search, the traversal search method is more computationally expensive than other DPD methods. Obviously, Weiss's method has a great improvement in the running time. Compared with Weiss's method, the run time of our method is further reduced (approximately 40%). However, the reduction of computational complexity of our method does not deteriorate localization performance significantly (also see Figure 3). As indicated by these finding, the proposed approach receives acceptable localization performance with a low computation resource cost. Thirdly, the CDF curves of all methods at the nominal constellation are depicted in Figure 4a,b at the SNR level of 5 dB and´10 dB, respectively. It can be readily observed that the traversal search curve is highest at a different localization error level. The proposed method curve is associated with Weiss's DPD curve, which indicates it can receive high localization accuracy with high SNRs. From Figure 4, the two-step method curve deviates from the DPD methods. The proposed method curve is slightly apart from the curve of Weiss's method, which is acceptable in terms of low complexity.  Finally, to further examine the influence of system parameter values on localization performance, we set SNR =´5 dB and exhibit the RMSE curves for K,N,L, and v in Figure 5, respectively. Additionally, the simulation conditions are K ranging from 2 to 10, N ranging from 50 to 200, L ranging from 1 to 6, and v ranging from 100 to 1000. It can be readily found in Figure 5a-c that all algorithms perform better in terms of localization accuracy, as more system parameter information exists. These results imply that the increase of these system parameters can enhance the localization performance. However, more parameter requirements mean greater system overhead and computational pressure. Consequently, our method can play an excellent role in the presence of limited computational resources. Furthermore, Figure 5d indicates that v does not have much impact on position accuracy.

Conclusions
The traditional two-step methods have high computational efficiency but low localization performance. To enhance positioning accuracy, DPD approaches have been developed. Since the DPD technology locates the target from the signal directly, computational complexity of this estimator is high. To solve this problem, this paper proposes a fast ML-based direct localization method using an EM algorithm based on time delay and Doppler shift. The EM scheme was developed to solve the ML problem in the DPD model. We simplify this ML problem via the theories of Gibbs inequality and Laplace approximation, transforming the high-dimensional nonlinear search into a few one-dimensional searches. Simulation results show that the proposed algorithm has operates faster than other DPD methods. Furthermore, when limited computing resources appear, our method leads a more efficient result on the basis of guaranteeing high localization accuracy. Therefore, compared with other localization approaches, the proposed method becomes a good way to balance localization accuracy and computational complexity.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Evaluation ofx l,k andḠ l,k via Laplace Method
We consider the multivariate functions h pzq : R n Ñ R m and f pzq : R n Ñ R. An integral is expressed by Assume that f pzq has a global maximum within the interval. According to [27], the Laplace approximation is written as where F pz˚q is the Hessian of f pzq at z˚. Next we approximate the integral ofx l,k andḠ l,k . For the sake of brevity, we omit l, k in the subsequent derivation. The expression of the conditional expectation is The polar coordinates can be written as Similarly, we haveḠ « G pr˚q.

Appendix B. Derivation of the Cramér-Rao Bound
It is well known that the CRB provides a benchmark for the best localization accuracy for any unbiased estimator. This section derives the compact CRB formula for source position with known signal waveforms and unknown transmitted time.
Let D l,k pθq " b l,k C l,ksk . According to [22], the q, jth entry of the fisher information matrix of unknown parameter θ is shown as rJs q,j " (A10) For easy derivation, θ can be redefined by Since submatrices of Equation (A12) with respect to diagonal symmetry are transposed matrices of each other, we can easily obtain the remaining submatrices.
Thus, the inverse of the CRB of position can be obtained bỳ CRB p˘-1 =J pp -" J pb pRq J pb pIq J pt 0 ı » -- This concludes the derivation.