Direct Position Determination of Unknown Signals in the Presence of Multipath Propagation

A novel geolocation architecture, termed “Multiple Transponders and Multiple Receivers for Multiple Emitters Positioning System (MTRE)” is proposed in this paper. Existing Direct Position Determination (DPD) methods take advantage of a rather simple channel assumption (line of sight channels with complex path attenuations) and a simplified MUltiple SIgnal Classification (MUSIC) algorithm cost function to avoid the high dimension searching. We point out that the simplified assumption and cost function reduce the positioning accuracy because of the singularity of the array manifold in a multi-path environment. We present a DPD model for unknown signals in the presence of Multi-path Propagation (MP-DPD) in this paper. MP-DPD adds non-negative real path attenuation constraints to avoid the mistake caused by the singularity of the array manifold. The Multi-path Propagation MUSIC (MP-MUSIC) method and the Active Set Algorithm (ASA) are designed to reduce the dimension of searching. A Multi-path Propagation Maximum Likelihood (MP-ML) method is proposed in addition to overcome the limitation of MP-MUSIC in the sense of a time-sensitive application. An iterative algorithm and an approach of initial value setting are given to make the MP-ML time consumption acceptable. Numerical results validate the performances improvement of MP-MUSIC and MP-ML. A closed form of the Cramér–Rao Lower Bound (CRLB) is derived as a benchmark to evaluate the performances of MP-MUSIC and MP-ML.


Introduction
We are interested in high precision positioning for shortwave signal sources in this paper. Two-step methods, such as the Angle Of Arrival (AOA) method, were usually used for shortwave signal positioning, and the methods provided a poor performance in a low Signal Noise Ratio (SNR) scenario. It has been shown that available prior knowledge on deterministic multi-path components can be beneficial for localization [1]. Kietlinski-Zaleski Jan presented techniques to benefit from signal reflections from known indoor features such as walls [2]. Inspired by those ideas, we propose a novel geolocation system architecture to locate the shortwave sources. This new architecture, termed "Multiple Transponders and Multiple Receivers for Multiple Emitters Positioning System (MTRE)", uses multiple transponders and receivers with known locations to locate multiple narrow band signals. The raw signals are transferred "in band" (i.e., as a man-made multi-path) by the transponders, and there is no need for a network infrastructure or an out-of-band channel bandwidth, which are required in an up/down converter system. In order to avoid the interference between the receiving and the sending signals of a transponder, we use different polarization modes to isolate the signals. In an MTRE system, man-made multiple paths from an emitter to a receiver are made to improve the positioning precisions and extending the positioning range.
Multi-path propagation is a really important problem in outdoor and indoor positioning systems, and it is still the main source of estimation errors for range-based indoor localization approaches [3,4]. The recent research in dealing with multi-path either tries to detect these situations statistically based on the received signals [5,6] or to directly mitigate the corresponding errors with statistical techniques [7,8]. Some algorithms for indoor localization make use of, e.g., the cooperation of multiple agents to overcome multi-path situations [9]. Arrays were used for beam-forming to separate signals from different directions, and the multi-path positioning problem is simplified into a single path positioning problem [10,11]. Furthermore, location fingerprinting, e.g., Received Signal Strength (RSS)-based methods, has been widely used in harsh environments [12,13]. It makes use of a priori training signals in multiple regions of the environment to train a classification algorithm [14]. However, the required training phase, as well as the missing flexibility w.r.t. changes in the environment may limit its application.
Most of the above literature focused on the UWB signals and indoor positioning applications, and two-step approaches were adopted to locate the emitters. The very high bandwidth of the UWB signal translates into very good time resolution and makes the UWB signal resistant to multi-path. It is possible to extract parameters, e.g., RSS, Time Of Arrival (TOA), AOA, Time Difference Of Arrival (TDOA) and Frequency Difference Of Arrival (FDOA), from UWB signals in the presence of multi-path propagation and to locate the emitter based on those parameters [15,16]. However, narrow-band systems have a low time resolution, and it is difficult to get the measurements in the first step.
The Direct Position Determination (DPD) methods were proposed in [17] for single narrow band signal positioning and in [18] for multiple narrow band signal positioning. A DPD approach collects data at all sensors together and uses both the array responses and the Times Of Arrival (TOA) at each array, in contradiction with the two separate steps: parameter measuring and location determination. From the optimization theory point of view, two-step methods are sub-optimal, since the parameter estimation in the first phase is done independently, without considering the constraints that the measurements must correspond to the same source position. DPD methods overcame the problem of associating estimated parameters with their relevant sources and was shown to outperform two-step methods, especially in low SNR scenarios [19].
There have been only a few attempts to improve the accuracy of emitter positioning in the presence of multi-path propagation under the DPD framework. Most of the existing DPD methods were developed for a single-path channel in which the multi-path was modeled as additive noise [20]. In [20], the single path DPD was tested with a channel with two paths scenario and showed improved performance over two-step methods. The DPD with small local scattering was studied in [21,22]. In a scattering scenario, sensors were affected by a set of virtual emitters, which were placed randomly in close proximity to the real emitter. It was assumed that the positions of the virtual emitters were i.i.d, and each position forms a 3D Gaussian distribution. Odel Bialer and Dan Paphaeli and Anthony J. Weiss [23,24] proposed a positioning algorithm for a dense multipath environment. Each received signal was obtained by convolving the transmitted pulse with a channel impulse response, and only the first arrival cluster (the direct path) was taken into consideration in their work. The signals reflected from other objects were not modeled in their work.
Papakonstantinou and Slock [25,26] considered a simplified single-bounce multipath model. The model assumed that the transmitted signal did not bounce over more than one scatter. They jointly estimated the position of the target and scatters. They studied the single emitter positioning problem in the presence of multi-path propagation and assumed that the waveform of signal and path attenuations were known in advance. Miljko and Vucic [27] proposed a novel direct geolocation of an Ultra WideBand (UWB) source in the presence of multi-path using the MUltiple SIgnal Classification (MUSIC) method and focusing matrices. Only one emitter was taken into consideration, and the path attenuations are known in advance in their work. Bar-Shalom et al. [28,29] proposed a transponder-added Single Platform Geolocation (SPG) model. A single emitter and single receiver were assumed in the SPG model. They stated that the SPG model achieved a similar performance to the multiple-RX DPD algorithm. The multiple-RX DPD algorithm mentioned in their works assumed that the transponders were replaced by the receivers directly. In a weak signal location application, a single receiver cannot receive signals from all transponders stably. Some paths may be blocked or disrupted. Multiple emitters, multiple transponders and multiple receivers need to be taken into consideration in a weak signal positioning application.
All unknown parameters should be estimated together in a DPD model, and this leads to a large-scale parameter searching. MUSIC methods calculate the spatial spectrum of each candidate position rather than the combinations of all emitter positions. Amar et al. [30] studied that multiple known and unknown radio-frequency signals under the LoS (Line of Sight) channel assumption. A simplified MUSIC algorithm was adopted to avoid the large-scale parameter searching. The cost function in [30] maximized the projection of the array manifold onto the signal subspace rather than minimizing the projection onto the noise subspace. The simplified cost function took advantage of the maximization of the convex Quadratic Programming (QP) with linear constraints, and the eigenvalue structure was adopted to avoid the searching of path attenuation parameters in their work. The simplified MUSIC worked well in an LoS propagation context, but it had a poor performance in a multi-path propagation scenario due to the singularity of the array manifold. Minimizing the projection of the array manifold onto the noise subspace overcomes the shortcomings of a signal subspace projection method. However, the eigenvalue system fails to resolve the minimization programming.
Existing DPD methods mainly focus on narrow-band signal positioning [18,29,[31][32][33] and usually assume that the carrier phase does not carry the propagation delay information. Complex channel attenuations were estimated to eliminate the influence of carrier phase misalignment in a narrow-band signal positioning method. We point out that the narrow-band assumption will lose the phase information in an LoS positioning application, and it will not be able to locate emitters at all in a multi-path positioning application. We add constraints that path attenuations are nonnegative real numbers in our model. In an existing DPD model, path attenuations are complex numbers and have only one equation constraint (the norm of path attenuations is one), and the Lagrange-multiplier method is very effective at solving the optimization with equation constraints [34]. However, it is difficult to solve an optimization with inequality constraints (the path attenuations should be greater than zero). We are required to design an efficient algorithm to solve the QP with inequality constraints.
The performance of a MUSIC method is determined by the precision of the covariance matrix estimation. In a time-sensitive application, the number of snapshots is not enough, and it is difficult to estimate the covariance matrix precisely. The maximum likelihood method maximizes the likelihood function of the received data rather than estimating the covariance matrix, and it achieves a better performance than that of the MUSIC method. However, the dimension of the searching space turns out to be unacceptable in the maximum likelihood method.
Our motivation is to develop a simple and accurate positioning model and corresponding algorithms for the case of unknown waveform signals and multi-path environment. We establish a Multi-path Propagation (MP)-DPD model for the scenario of multiple emitters, multiple transponders and multiple receiving arrays. It can be viewed as a modified and extended version of the SPG model proposed in [29]. The MP-DPD reduces the risk of paths being blocked or disrupted and fixes the constraints on path attenuations. Multiple emitters can be simultaneously positioned in the MP-DPD model, as well. MP-MUSIC and MP-ML methods are proposed to reduce the time consumption of the optimization. The numerical results and the Cramér-Rao Lower Bound (CRLB) analysis show that the MP-MUSIC method has a lower computing complexity than MP-ML, especially in the case of a complex multipath scenario. The MP-ML method is more precise than MP-MUSIC, especially in the case of positioning with limited snapshots. An Active Set Algorithm (ASA) for the MP-MUSIC and an iterative algorithm for the MP-ML are developed to reduce the computational complexity of the methods further. Numerical results demonstrate that the MP-MUSIC and MP-ML proposed in this paper outperform the conventional methods.
The paper is organized as follows: Section 2 outlines the problem formulation, and an MP-DPD model is established in this section. The MP-MUSIC method, the MP-ML method and corresponding algorithms are proposed in Sections 3 and 4, Numerical performance examples of these algorithms are given in Section 5. The final conclusions are given in Section 6. Finally, the detailed descriptions of the ASA algorithm, the iterative algorithm for MP-ML method and the derivation of the CRLB are provided in the Appendix.

Problem Formulation
Consider that there are D emitters located at p e = [p T e (1), p T e (2), ..., p T e (D)] T and L passive transponders placed at The signals transmitted by the emitters are reflected by the transponders and intercepted by N receiving arrays. Each array includes M antennas. The centers of the arrays are located at It is assumed that the locations of the transponders and the receiving arrays are known a priori and that the signal waveforms are unknown. The scenario is depicted in Figure 1.  Denote the signal propagation delay between the d-th emitter and the -th transponder byτ d . Denote:τ whereτ nm is the propagation delay between the -th transponder and the m-th antenna in the n-th receiving array.τ n is an M × 1 column vector, which represents the propagation delays from the -th transponder to the n-th receiving array.τ n is known in advance, and it is independent of the emitter positions.
The path attenuation from the d-th emitter to the n-th receiving array, which is reflected by the -th transponder, is denoted by α d n . The path attenuation coefficients are assumed as non-negative real numbers, and the rationality of the assumption will be discussed in detail in Section 3.1.2. We assume that the antennas in a receiving array are uniform, and all antennas in an array share the same path attenuation coefficient.
The time-domain model of the signals that are received by the n-th receiving array is: wherer n (t) is an M × 1 column vector, which represents M snapshots at time t of the n-th receiving array.s d (t) is an M × 1 column vector, which represents M snapshots of the d-th source signal at time vector t t −τ n −τ d − t d .n(t) is an M × 1 noise vector at time t. 0 ≤ t ≤ T, and t d is the unknown transmit time of the emitter d. We assume that the path attenuation, α d n , remains constant during the observation time interval. This paper mainly focuses on the positioning of deterministic, but unknown signals. It is assumed that source signals are independent of one another, and there is no further requirement for the code or waveform of the signals. The frequency-domain model for the k-th DFT coefficients is given by: where:ã where s d (k) is the k-th Fourier coefficient of the d-th source signals d (t), t ∈ [0, T]. r n (k) and n(k) are M × 1 vectors of the k-th Fourier coefficients ofr n (t) andn(t).ã n (k) is an M × 1 vector, which denotes the generalized array response of the n-th receiver at frequency ω k . Make (3) into matrix form: where: where ⊗ is the Kronecker product and I N is an identify matrix with a size of N × N.
Denote the second moments of variables by: E{n(k)n(k) T } = 0, where I MN is an identity matrix with a size of MN × MN. R(k) is a covariance matrix of received signals at frequency ω k . σ is the noise standard deviation. The observed signal of each antennar(t) is partitioned into J sections, and each section is Fourier transformed. The k-th Fourier coefficient of the j-th section is denoted by r j (k). The covariance matrix at frequency ω k is estimated by: The relationship between the received signals and emitter positions has been established in (5), and it is named the MP-DPD model. The MP-DPD model optimizes the emitter positions directly to achieve a more accurate estimation. Two DPD methods are proposed under the MP-DPD framework: The array manifold projection onto the noise subspace is adopted as the cost function in the MP-MUSIC method, and the likelihood function of the received signals is adopted in the MP-ML method. If the number of snapshots is sufficient, MP-MUSIC consumes less time than MP-ML without degrading the performance. Besides, if the number of snapshots is not enough, MP-ML obtains more precise position estimations than MP-MUSIC.

MP-MUSIC Method
We analyze the shortcomings of the existing MUSIC method in the multi-path propagation positioning firstly and establish an MP-MUSIC model that is suitable for the multi-path environment positioning. Finally, the corresponding algorithm is given at the end of this section.

The Limitation of Existing MUSIC Methods
We introduce the Signal Subspace Projection MUSIC (SSP-MUSIC) method, which was commonly used in the DPD model firstly, and develop a Noise Subspace Projection MUSIC (NSP-MUSIC) method to overcome the shortage of SSP-MUSIC. Finally, we discuss the performances of SSP-MUSIC and NSP-MUSIC, which were adopted in a multi-path positioning application.

SSP-MUSIC
Alan Amar and Anthony J. Weiss studied the positioning problem of multiple unknown radio-frequency signals in [30]. The MUSIC method for the LoS propagation positioning was proposed in their works. Alan Amar maximized the manifold projection onto the signal subspace rather than minimizing the projection onto the noise subspace. The programming model of Amar's method was defined as: [p e ,α] = arg max F(p e , α) = α H D(p e )α, where: where α is the path attenuation vector, C ND is the set of complex column vectors with the length of ND, · F is the Frobenius norm of a matrix, p e is the vector of emitter positions, I N stands for the N × N identity matrix, 1 M stands for an M × 1 column vector of ones, M stands for the number of antennas of each array, Γ(k) is the array manifold matrix at frequency ω k , N is the number of receivers, D is the number of emitters, K is the frequency points of received signals and U s (k) is made up of the eigenvectors of the covariance matrix of received signals corresponding to the D largest eigenvalues. The other parameter notations can be found in [30]. Since (10) is a quadratic convex optimization with linear constraints in the complex field, the maximum of the cost function is the maximal eigenvalue of the matrix D(p e ) [31]; thus, the optimal cost function value is F * (p e ) = λ max {D(p e )}, where λ max {·} represents the maximal eigenvalue of a matrix. Benefiting from the simplified cost function and the eigenvalue system, Amar's method reduced the searching dimension from 2DK + 2(N − 1)D + 3 to three.

NSP-MUSIC
A noise subspace projection MUSIC method is proposed in this paper to remove the simplification of the SSP-MUSIC. NSP-MUSIC minimizes the manifold projection onto the noise subspace: Reorganize the items in (14): [p e ,α] = arg maxF(p e , α) = 1 where: and D(p e ) has been defined in (12). SSP-MUSIC is viewed as a simplified version of (15). In a direction finding application, it is assumed that the path attenuation of each antenna in an array has been normalized in advance, and α H I(p)α is a known constant item. The dropping of the constant item in an optimization is reasonable. However, in a multi-path positioning application, α H I(p)α changes with the change of α, and the dropping of α H I(p)α is unreasonable. Following the standard noise subspace MUSIC method and the model for the multi-path positioning, we define the cost function of NSP-MUSIC: where I MN is an identity matrix with a size of MN × MN. U s (k) is a matrix consisting of the eigenvectors of R k corresponding to the D largest eigenvalues. p e and α are decision making variable vectors representing the candidate emitter positions and the corresponding path attenuations. In order to facilitate a unique solution, we assume that the norm of α is one. a(k) is the array manifold vector for p e at frequency ω k . Unfortunately, the cost function requires an LN − 1 + 3 dimensional searching, and it is difficult to get the optimal solution over such a high dimensional space. Note that the d-th column of matrix A(k) in (5) is denoted by: where α d [α T d1 , α T d2 , ..., α T dN ] T represent attenuations of paths from the d-th emitter. The vector a(k) of a candidate emitter in the MUSIC algorithm (17) is similar to a d (k), but the position of the d-th emitter is replaced by the candidate emitter position p e . Denote Γ(k) Γ d (k) to simplify the explanation, and substitute (18) into (17): Tom Tirer and Anthony J. Weiss studied similar programming in [35]. They transformed the cost function into: where λ min {·} represents the minimal eigenvalue of a matrix. It is a promotion result of the maximization QP, but it is only a "not bad" solution rather than the optimal one. If Γ(k) is singular, E(p e ) turns out to be singular, and λ min {E(p e )} = 0. In this case, the cost function reaches a peak. NSP-MUSIC only finds the solutions that make Γ(k) singular rather than true emitter positions. From another point of view, if ∃i, j, which satisfy e i (p e ) ≈ e j (p e ), where e i (p e ) and e j (p e ) are the column i and column j in E(p e ), the matrix E(p e ) turns out to be singular or near singular (It should be noticed that, in a single path positioning application, Γ(k) is a block diagonal matrix, and each block is an M × 1 column vector. It is impossible that Γ(k) has the same two columns, but this is possible for a multi-path positioning application. In a multi-path environment, Γ(k) is a block diagonal matrix, and each block is an M × L matrix. It is possible that the block is singular.). In this case, the optimal estimations of path attenuations areα = [α 1 ,α 2 , ...,α ·n ] T , where: α z is a feasible solution that satisfies (20). Substitute (22) into (19); Q(p e , α) → +∞. If α are complex scaled path attenuations, the cost function of NSP-MUSIC will reach a peak where E(p e ) is singular or near singular.
Above all, if the manifold matrix Γ(k) is singular or near singular, SSP-MUSIC will fail to get the emitter positions. Besides, if Γ(k) is singular or near singular and α are complex path attenuations, NSP-MUSIC also fails to get the emitter positions. In the next section, we will discuss the singularity of the manifold matrix Γ(k) and the necessity for non-negative real number constraints for path attenuations.

Singularity of the Manifold Matrix in the Presence of Multi-Path Propagation
We have discussed that SSP-MUSIC and NSP-MUSIC will fail to locate the emitters if Γ(k) is a near singular matrix, and we get the conditions of a candidate that makes Γ(k) near singular in Appendix A.
Unfortunately, Γ(k) will always be near singular. For example, in a shortwave positioning application, the size of a shortwave antenna is large. If the receivers need to be installed on mobile platforms (e.g., aircraft), only one antenna can be installed in a receiver (M = 1). In this case, from the condition in Theorem A2 in Appendix A, a 1 n (k) = a 2 n (k), k = 1, 2, ..., K always are satisfied, and the manifold matrix Γ(k) of a CSMCis a singular matrix.
In another application, the transponders are installed on a mobile platform (e.g., Unmanned Aerial Vehicle (UAV) platform or satellite platform). If one transponder is relatively close to another transponder, or two transponders and a receiving station are near collinear, this makes a 1 n (k) ≈ a 2 n (k), k = 1, 2, ..., K. In this situation, the conditions in Theorem A2 in Appendix A are satisfied, and Γ(k) turns out to be near singular.
In addition, it is necessary to down convert the Radio Frequency (RF) signals to baseband signals firstly to avoid the multi-peak searching of the cost function (see Figure 2). The suboptimal peaks in Figure 2a are caused by the carrier wave. The higher the frequency of the carrier, the more suboptimal the peaks in the cost function. Besides, there is only one peak in the cost function for a baseband signal positioning model (see Figure 2b). Figure 2b is an up envelopeof Figure 2a. If the cost function is a surface with multiple peaks, it is difficult to develop a searching strategy except for a high density grid searching. However, it is easy to get the global optimal solution for a continuous function with a single peak (e.g., Steepest Descent Method (SDM), Newton Method (NM)). In a baseband signal positioning application, λ R, where λ is the wave length corresponding to the maximal frequency of the baseband signal, and R is the radius of the circle receiving array. If λ R, it is easy to satisfy a 1 n (k) ≈ a 2 n (k), k = 1, 2, ..., K, and Γ(k) is near singular.

Non-Negative Real Path Attenuation Constraints
Our model requires that path attenuations must be non-negative real numbers, but complex values in existing studies for narrow band signal positioning [18,29,[31][32][33]. Weiss ignored the path attenuations (set the attenuation α = 1) in wide-band emitter positioning in [36].
In a narrow band signal positioning application, it is assumed that a n (k) ≈ a n (k 0 ) a n , where k = 1, 2, . . . , K, and ω k 0 is the carrier frequency. Based on the above assumptions, (3) turns out to be: where e jω k 0 τ is a phase adjustment item to satisfy e jω k 0 τ a n (k 0 ) → a n (k). Existing studies used the envelope information only to estimate the propagation delay and dropped the carrier phase information. e jω k 0 τ was an adjustment factor for carrier phase alignment, and it was used to reduce the interference with the propagation delay estimation. Denoteᾱ d n α d n e jω k 0 τ ; (23) turns out to be: whereᾱ d n is a complex scalar representing the "channel attenuation" (It is not a real channel attenuation coefficient, but an equivalent parameter, which is determined by the real path attenuation and the model error caused by the narrow band signal assumption), and a n denotes the generalized array response matrix. The commonly-used DPD methods [18] modeled the received signal as: (25) is an LoS positioning model for narrow band signal positioning, while (24) is an NLoS positioning model. Existing models with complex path attenuation assumptions are viewed as simplified models of (3) in a narrow band signal positioning application. We point out that the simplified model (24) can not be adopted either in a MUSIC method or in an ML method in a multi-path propagation and unknown wave form application. We have obtained that the manifold matrix Γ(k) may be singular in Section 3.1.3 and discussed that if the path attenuations were complex numbers, a MUSIC method could not obtain the emitter positions correctly in Section 3.1.2. We will discuss this further in Section 4.1 to explain the necessity of the real and non-negative constraints in an ML method. Overall, we develop (3) as the signal model for positioning and constrain the path attenuations in R + .

Mathematical Model of MP-MUSIC
In an MP-MUSIC method, the optimal estimation of α for a fixed emitter position p e is given by solving the following programming: where · 1 is the one norm of the vector (that is, the sum of the absolute values of the elements of the vector) and I MN is an MN × MN identity matrix. The programming is a non-linear programming with real value constraints. There is an LN dimensional searching for a candidate emitter position, and it is difficult to solve the programming directly. We remove the imaginary items in the programming without changing the optimal solution firstly and prove the convexity of the modified programming later. An iterative algorithm named ASA is proposed to solve the convex programming after the proof.

Remove the Imaginary Items in the Programming
Rewrite the objective function by: where Ψ Re[E(p e )] is the real part of E(p e ) and Φ Im[E(p e )] is the imaginary part of E(p e ). Φ is a Hermitian matrix, which satisfies: where Φ i,j is the element at the i-th row and j-th column of the matrix Φ. Since α is a non-negative real value vector and the diagonal elements of Φ are zeros, Substitute (29) into (27), and the cost function turns out to be: The programming (30) is a QP in the real field. If there are no inequality constraints and the objective function is convex, the Lagrange multiplier method is effective for solving the programming. However, the non-negative constraints of path attenuations are necessary due to the singularity of the array manifold. To obtain the optimal solution of (30), we verify the convexity of the programming firstly and design an algorithm to solve the convex programming.

Theorem 1. (30) is a convex quadratic programming with linear equality constraints and lower bounds.
Proof. where: where U n is the noise subspace of the received signal. ∀x ∈ R LN , Since ∀x ∈ R LN , x H Φx = 0, and substituting (33) into (34), we get x H Ψx ≥ 0, and the objective function is Positive Semi-Definite (PSD). The optimization problem (30) is a convex quadratic programming with linear equality constraints and lower bounds.
It is possible to find the global optimal solution for a convex quadratic programming [37,38]. The interior-point algorithm or any Heuristic Searching Algorithm (HSA) can be adopted to solve the optimization problem with equality constraints and lower bounds. However, those algorithms apply numerical searching strategies with low efficiencies. We introduce a faster algorithm named the Active Set Algorithm (ASA) in this paper to obtain the global optimal solution base on some theoretical analysis.

Active Set Algorithm
The first widely-used algorithm for solving a similar problem is the active set method published by Lawson and Hanson [34,39]. They proposed an active set method to solve the Non-Negative Least Squares (NNLS).

Definition 1.
Active set [40] In mathematical optimization, a problem is defined using an objective function to minimize or maximize and a set of constraints: Given a point x in the feasible region, a constraint: is called active at x if g i (x) = 0 and inactive at x if g i (x) > 0. The set of active ones is called the active set and denoted by A(x) = {i|g i (x) = 0}.
We describe the active set method for solving the quadratic programs of the form (30) containing equality and inequality constraints based on the methods described in [34,39].
Denote the optimal solution of (30) byα. If the active set of the optimal solution A(α) were known in advance, we could find the optimal solutionα by applying techniques, such as the Lagrange multiplier method, for equality-constrained QP. The prior knowledge of the active set accelerates the algorithm effectively.
The searching path of the active set algorithm is strictly in the feasible region. Choose a feasible solution as the initial point of the algorithm. Solve the QP with equality constraints in the work set, and get the optimal searching direction. if the searching direction is blocked by some constraints not in the work set, add the constraints, which block the searching path firstly, into the work set. Resolve the new QP with the updated work set, until the searching direction is not blocked by any constraints. The algorithm reaches a local optimum for the current work set. To get an even better solution, we drop one active constraint in the work set to relax the programming. If the objective function cannot be decreased for all constraints in the work set, we get the optimal solution of the original programming. Otherwise, drop the constraint that causes the fastest decrease.
The detail of the ASA is described in Section 1 of the Supplementary File. The spatial spectrum of an emitter is determined by substituting theα into the MUSIC cost function (19).

MP-MUSIC Algorithm
The spatial spectrum of the emitter positions requires only a three-dimensional searching, and the size of E(p e ) is LN × LN, which is usually rather small. The detailed procedure of the MP-MUSIC algorithm is represented in Algorithm 1. The performance of the MP-MUSIC algorithm is determined by the estimation precision ofR(k). If the number of snapshots is not enough, neither the covariance matrixR(k) nor the spatial spectrum q(p e ) can be estimated precisely.
In a time-sensitive positioning application, it is difficult to get enough snapshots to estimateR(k). We develop a Maximum Likelihood method in the presence of Multi-path Propagation (MP-ML) to estimate the emitter positions directly.

MP-ML Method
The MP-ML method maximizes the conditional likelihood function of the received signals. The noise is assumed as Additive White Gaussian Noise (AWGN) with a known standard deviation σ,

Mathematical Model of MP-ML
The likelihood function of the received signals is: where r is the observed data, Σ is the covariance matrix of noises, which is defined in (6), and the unknown parameter vector θ [p T e , α T , s T ] T consists of: where α dn andš(k) have been defined in (5). The log-likelihood function of (37) is: Remove the constant items, and get the modified cost function of MP-ML: The searching space dimension of (40) is 3D + DNL + DK, and it is necessary to reduce the searching space dimension. For fixed attenuations α and emitter position combination p e in (40), the optimal estimation of the source signals at frequency ω k is: where (40), where η [p T e ,ᾱ T ] T ,ᾱ = αI D [α 1,1,1 , . . . , α d, ,n , . . . , α DLN ] T , I D is a column vector of D ones. P A (k) = A(k)A + (k) is the projection matrix of A(k). Expand (42): and move the constant items r H (k)r(k). Applying the properties of the projection matrix P A (k) = P H A (k) and P H A (k)P A (k) = P A (k), we get the modified programming of MP-ML: There are two differences between our model and Bar-Shalom Ofer's model in [28]. The first one is that only a single emitter and a single receiver were modeled in their work, but multiple emitters and multiple receivers are taken into consideration in our work. The second one is that there are complex path attenuations in [28], but real non-negative path attenuations in our model.
The cost function in [28] was modeled as: where f(k) Γ H (k, p e )r(k), C(k) Γ H (k, p e )Γ(k, p e ). Section 3.1.1 has discussed that Γ(k, p e ) may be singular, and C(k) may be singular, as well. When C(k) is singular, there is an α that satisfies α H C(k)α = 0, and the cost function Q(η) reaches the peak. However, the candidate emitter position is not the true emitter position. If Γ(k, p e ) is near singular, f(k) = Γ H (k, p e )r(k) is near singular, as well. The numerator and denominator of the cost function Q(η) both tend to zero, and the value of the cost function turns out to be unstable. The noise level will seriously affect the value of the cost function in this case, and the model cannot find the emitter accurately.
The searching space dimension of (44) has been reduced to 3D + DNL, but it is still difficult to solve such a high dimensional non-linear programming. We propose an iterative algorithm in this paper to get the estimation of path attenuations to reduce the time consumption of the MP-ML.

Remove Imaginary Items in the Programming
Substitute the constraints into the objective function of (44), where f(k) Γ H (k, p e )r(k), C(k) Γ H (k, p e )Γ(k, p e ). Henk A. L. Kiers studied a similar convex optimization problem in [41]. Ofer Bar-Shalom and Anthony J. Weiss study the complexity form of the optimization in [28] and its application in [29]. The programming in our work has the complex f(k) and C(k), but the decision making variables α are real non-negative values. We modify the iterative process in [28,41] to satisfy the real non-negative constraints in our work.
The cost function (46) can be rewritten by: where tr(·) is the trace operator of a matrix. Since C(k) and f(k)f H (k) are Hermitian metrics, ∀α satisfy: The complex matrices f(k) and C(k) are replaced by the matricesf(k) andC(k) with real number elements: The complex non-linear programming with real constraints (46) is simplified to be a real non-linear programming (49).

An Iterative Algorithm for Solving MP-ML
We introduce a theorem firstly and then give an iterative algorithm for solving the programming (49). Theorem 2.ᾱ i is a feasible solution of (49), and a better solution of (49) is obtained by solving the following programming:ᾱ s.t.ᾱ ≥ 0. where: I L is an identify matrix with a size of L × L and I N is an identify matrix with a size of N × N.Ū is the Singular Value Decomposition (SVD) of the following item: The proof of Theorem 2 is given in Appendix B. The programming (50) is a linear least squares with bound constraints, and the Trust-Region-Reflective (TRR) algorithm is adopted to solve the programming. The detail of TRR is described in [42][43][44].
Following Theorem 2 and (A11) in the Appendix B: We get an iterative method to obtain a better solution than the previous step. Denote the initial estimation of the unknown channel attenuations by α 0 . The solution of (50) gives a better estimation of the unknown parameter vector due to the inequality of (51). Furthermore, we get a better estimation α 2 , so that H(ᾱ 2 ) ≤ G(α 1 ,ᾱ 2 , p e ) ≤ G(α 1 ,ᾱ 1 , p e ) = H(ᾱ 1 ). Thus, H(ᾱ 2 ) ≤ H(ᾱ 1 ) ≤ H(ᾱ 0 ). The detail of the iterative procedure is shown in the Algorithm 2 in the supplementary file.

Getting the Initial Value
The performance of the iterative algorithm is determined by the initial valueᾱ 0 . The path attenuations from the MP-MUSIC algorithm are used as the initial value of the MP-ML.
The initial path attenuations of the emitter d are denoted byᾱ 0 (d), and they are estimated by the MP-MUSIC in Algorithm 1.
For a fixed emitter position combination p e , evaluate the MP-MUSIC algorithm to get the initial path attenuationsα d of the position p e (d), d = 1, 2, ..., D.
Reshapeα = [α T 1 ,α T 2 , , ...,α T D ] T to get the initial attenuations vector: where:α andα dn is estimated from MP-MUSIC. while p e ∈ p do Evaluate: F(p e ) ← MP_MUSIC(r,p e ); p e ← p e + ∆ p ; end Set: q * = +∞ i = 0; Get the initial p 0 e from the spatial spectrum F(p e ), p e ∈ p; return q(p e ) end

Numerical Examples
Some numerical examples are given to demonstrate the performances of the above algorithms.  The path attenuation coefficients are drawn from a uniform distribution between zero and one. The SNR is defined in terms of "post-processing SNR", which is given by:

Scenario Setting and Performance Index Definition
Root-Mean-Squared Error (RMSE) of the estimated position is adopted as the performance index of the algorithms. RMSE is given by: where N s is the number of Monte Carlo runs, D is the number of emitters, p e are the real emitter positions andp e are the estimated emitter positions A scalar quantity of the CRLB matrix corresponding to RMSE is defined as: where λ max [CRLB(p e )] is the maximal eigenvalue of CRLB(p e ). The computation of the CRLB matrix is presented in Section 3 of the Supplementary File [45][46][47].

Performances of the MP-MUSIC Method
We compare the performances of the following MUSIC algorithms:

SSP-MUSIC and NSP-MUSIC in a Single Path Propagation Positioning Scenario
The single emitter is placed at [0,0] (km), and only the direct paths are taken into consideration. Figure 4 gives the spatial spectrum in a single path propagation positioning scenario (Actually, the spatial spectrum of a 3D positioning is a 3D spectrum, but in order to show the form of spatial spectrum more intuitively, the spatial spectrum displayed in this paper only gives a horizontal slice of the 3D spatial spectrum at the real Z value (where z = 0)).  We find the peak in the spatial spectrum of SSP-MUSIC (Figure 4a), but the NSP-MUSIC proposed in this paper (Figure 4b) holds a sharper peak and more precise estimation than SSP-MUSIC in a single path positioning context.
The following simulations are designed to study the performance of those MUSIC methods in the presence of multi-path propagation.

SSP-MUSIC, NSP-MUSIC and MP-MUSIC in a Multi-Path Propagation Positioning
We design four numerical simulations in this section. The simulation parameters are set as in Table 1, where R is the radius of the UCA, λ and f are the carrier wave length and frequency, M is the number of the antennas in a receiving array, B is the bandwidth of the source signals, K is the number of frequencies and J is the number of sections of a received signal. Baseband signal positioning: Figure 5 gives the spatial spectrum when R λ. The emitter positions p e that make E(p e ) singular or nearly singular constitute the yellow hyperbolic curves in the SSP-MUSIC spectrum and the NSP-MUSIC spectrum. Neither SSP-MUSIC nor NSP-MUSIC find the emitters correctly when R λ. We cannot find any peak in SSP-MUSIC. Three peaks are found in the spectrum of NSP-MUSIC, but they are disrupted by the hyperbolic curves. MP-MUSIC with real and non-negative constraints finds three sharp peaks correctly. Two transponders are close: Figure 6 gives the spatial spectrum when a transponder is close to the anther. The nearest distance of transponders is 20 km in the simulation. When two transponders are close, the array responses of the receiving array with respect to the two transponders are almost same. The yellow hyperbolic curves in the SSP-MUSIC and the NSP-MUSIC spectrum are the candidate positions, which make Γ(k) singular. SSP-MUSIC and NSP-MUSIC failed to locate the emitters in this context, but MP-MUSIC gets the emitter positions correctly.
Single antenna of each receiving station: Figure 7 gives the spatial spectrum when M = 1. The array responders are defined as a n (k) = 1, and an SMC will make Γ(k) singular.
General scenario: Figure 8 gives the spatial spectrum of a general parameter setting scenario. Although there is no deliberate construction of conditions that leads to Γ(k) singularity in the general scenario, SSP-MUSIC and NSP-MUSIC still cannot obtain the emitter positions. However, MP-MUSIC obtains the three emitter locations accurately.

Performances of MP-MUSIC-ASA and MP-MUSIC-IPA
We down convert the radio frequency signals to baseband signals firstly to avoid the multiple peak searching of the radio frequency signal positioning. The simulation parameters are set as in Table 2. The RMSE of MP-MUSIC-ASA, MP-MUSIC-IPA and the CRLB are given in Figure 9.
Since ASA could find the global optimal solution, the performance of MP-MUSIC-ASA should be no worse than that of MP-MUSIC-IPA. The numerical simulation results show that both MP-MUSIC-ASA and MP-MUSIC-IPA find the global optimal solution and have the same RMSE because of the convexity of the cost function.
Simulations are done on a server with Intel Xeon CPU E5-2630 v4, 16 G memory, and Matlab2016a. The MATLAB function "lsqlin" is adopted to verify the performance of IPA. We repeat the simulation 1000 times to get the distribution of the time-consumption of the ASA and IPA sections. It is assumed that the time consumptions t ASA and t IPA follow normal distributions, and t ASA ∼ N(5.8175 × 10 −5 , 4.8311 × 10 −6 ), t IPA ∼ N(3.4756 × 10 −3 , 2.2225 × 10 −5 ) (seconds).
Benefiting from the convex properties and a reasonable initial value of α in ASA, ASA consumes only 1.67% more time than IPA and has a more stable time consumption than IPA (It should be noticed that the time consumption of MP-MUSIC-ASA will be determined by the path attenuations and SNR. A low SNR and small α n make the constraint α n ≥ 0 an active constraint (α n = 0.). The time consumption of MP-MUSIC-ASA is deeply affected by the number of elements in the active set.

The Performance of MP-ML and MP-MUSIC
We compare the performances of MP-ML and MP-MUSIC in cases of different numbers of snapshots.

Insufficient Snapshots
If the number of snapshots is insufficient, we cannot obtain a reliable covariance matrix estimation of observations, and MP-MUSIC will fail to get the emitter positions. The simulation parameters are set as in Table 3. Table 3. Parameter setting in insufficient snapshot scenarios. The performances of MP-MUSIC and MP-ML with 32 snapshots (16 frequencies) are compared in Figure 10. Only 32 snapshots of each receiver are taken for positioning. Thirty two snapshots are not divided into sections to estimate the covariance matrix in NSP-MUSIC (that is K = 16, J = 1). The RMSE of NSP-MUSIC cannot be reduced further with the increasing of SNR because of the error of the covariance matrix estimation. The MP-ML obtains a better performance than the MP-MUSIC in the sense of insufficient snapshots.

Performances of Different K and J Combinations
We discuss the performances of MP-MUSIC with different K and J combinations in Figure 11. The total number of snapshots in the simulations is 256 (K · J = 128). The other parameters are set as in Table 4.   The performances of MP-MUSIC methods are effected by the error of covariance estimations and the bandwidth. A smaller J leads to a larger error of the covariance estimation; in addition, a smaller K leads to a smaller number of observation equations.

The Performances of Different Numbers of Snapshots
We compare the performances of the MP-MUSIC and the MP-ML with different numbers of snapshots in Figure 12 (The CRLB derived in the Appendix shows that the CRLB is determined by the signal snapshots. The CRLB in Figure 12 is the average of 100 simulations with random signal snapshots.). The other parameters are set as in Table 5.   Table 6, and the simulation results are given in Figure 13.   Figure 13b presents the performances corresponding to Figure 13a. When the snapshots are sufficient (K J = 128), the positioning accuracy improvement of MP-ML is not significant compared to MP-MUSIC, but it costs much more time to obtain the positions. In this case, applying MP-MUSIC to get emitter positions is a wise choice. When the snapshots are insufficient (JK = 16), MP-ML obtains a better performance than MP-MUSIC, although MP-ML costs much more time to get the results. It is worth much more time to obtain a much better performance when the snapshots are insufficient.
Benefiting from a good initial value of MP-MUSIC and an efficient searching strategy (steepest descent method), the time consumption of MP-ML is acceptable when the snapshots are insufficient. The time consumption of ASA in MP-MUSIC is determined by the a priori information of the initial active set W, and it can be derived with some other techniques, e.g., direction finding. MP-ML costs much more time than MP-MUSIC, especially in the the case of a large number of emitters. Fortunately, MP-MUSIC, which is adopted to obtain an initial value of MP-ML, provides "not bad" solutions quickly, and the iterative algorithm continues outputting better and better solutions (see Figure 14).
The simulation parameters are set the same as the second row in Table 6, and d = 3. MP-ML does not output the result until T 0 , because MP-MUSIC is running from zero to T 0 to obtain the initial solution of MP-ML. The RMSE of MP-ML is monotone decreasing and continuous for the output results. In real applications, we can find a trade-off between the time consumption and positioning accuracy to obtain an acceptable solution.

Performance of SGP and MP-ML
Single Platform Geolocation (SPG) mentioned in [28] used only one platform to position the single emitter; while MP-ML uses multiple receiving arrays and locates multiple emitters simultaneously.

SPG and MP-ML with a Single Emitter
Assume that there is only one emitter (D = 1, p e = [0, 0, 0] T km) in the RoI. Figure 15 compares the performance of SPG and that of MP-ML with different numbers of receivers. The other parameters are set as in Table 7.  We compare the performances of MP-ML with N receivers, N = 1, 2, 3, 4, and give the corresponding CRLB, as well. When D = 1 and N = 1, MP-ML degenerates into SPG. The RMSE of N receivers is obtained from the average of the RMSE of all combinations of N receivers in four positions, and the positions of the four receiving stations are defined as Layout A. The simulation demonstrates that the performance is significantly improved as the number of receiving stations increases. The program with four receiving stations gains nearly a 10 5 performance improvement compared to SPG.

Performance of Positioning Multiple Emitters
MP-ML has the ability of positioning multiple emitters synchronously. We analyze the performances of different numbers of emitters in this section.
We  Figure 16. The other parameters are set the same as in Table 7.   Figure 16 only gives the RMSEs and CRLBs when D = 1, 2, 3, 4, since the RMSE turns out to be unstable and the CRLB turns out to be +∞ when D = 5. This can be explained by the algebra principle that the number of unknown emitters cannot be greater than the number of transponders.
The numerical simulations and CRLB results demonstrate that the performance of MP-ML is influenced by the number of emitters. If D L, the number of emitters has only a slight effect on the performance of MP-ML, and if D = L, the performance will decline significantly. If D > L, MP-ML cannot find any emitter at all. In the case of time-sensitive positioning, the number of snapshots is not enough. The maximum likelihood estimation algorithm for Multi-path Propagation positioning (MP-ML) maximizes the likelihood function rather than calculating the covariance matrix of the observations to avoid the requirement of a large number of snapshots. We designed an iterative algorithm and proposed the strategy of choosing an initial solution to accelerate the solving of the programming. Numerical simulation results show that MP-ML can approach the Cramér-Rao Lower Bound (CRLB) relative to MP-MUSIC with the same data length, but MP-ML requires more computation time than the MP-MUSIC method.

Conclusions
Furthermore, we discussed the performances of MP-ML with different numbers of receiving arrays and emitters. SPG mentioned in [28] is viewed as a degenerate version of MP-ML (where D = 1 and N = 1). The numerical results shows that it is worthwhile to increase the number of receiving stations in the sense of a weak signal, although MP-ML increases the hardware costs, communication overhead and computational complexity.
We compared the performances and time consumptions of MP-MUSIC and MP-ML by numerical simulations. MP-ML obtains a more precise position estimation than MP-MUSIC, and MP-MUSIC consumes less time than MP-ML. In a specific positioning application, we choose the appropriate method according to the number of snapshots, the precision requirement and the calculation ability.
MP-ML has the ability of positioning multiple emitters synchronously. If the number of emitters is far less than the number of transponders, the number of emitters has only a slight influence on the positioning performance. If the number of emitters is equal to the number of transponders, the performance will decline significantly. If the number of emitters is more than the number of transponders, MP-ML cannot find any emitter at all.
An MTRE system requires more receiving arrays, more transponders and more computing resources compared to a Single Geolocation Platform (SGP) or a Direction Finding System (DFS). However, an MTRE system can locate multiple emitters synchronously and provides a higher positioning accuracy than SGP and DFS. It is suitable for some cost-insensitive applications, such as military and national security applications.
the center of the n-th receiving array andτ i , i = 1, 2, is the propagation delay from the candidate position to the i -th transponder.
The difference of an SMC and a CSMC is that all the antennas satisfy the equation for an SMC, but only the centers of the array satisfy the equation for a CSMC.
Theorem A1. If p e is a CSMC, there are at least two paths from p e to a receiving array, which satisfy: D 1,2 n (p e ) = zλ, z = 0, 1, 2, . . . , where n is the receiving array index, 1 and 2 are the indexes of two transponders in two paths,D 1,2 n p e − p t ( 1 ) F + p r (n) − p t ( 1 ) F − p e − p t ( 2 ) F − p r (n) − p t ( 2 ) F is the difference between the two path lengths, p r (n) is the center of the n-th receiving array, z is an integer, λ is the Least Common Multiple (LCM) of {λ 1 , λ 2 , . . . , λ K } and λ k is the wave length of frequency ω k .
where λ k = 2πc ω k is the wave length of frequency ω k , c = 3.0 × 10 8 m/s is the light speed constant, D 1,2 n p e − p t ( 1 ) F + p r (n) − p t ( 1 ) F − p e − p t ( 2 ) F − p r (n) − p t ( 2 ) F is the difference between the two path lengths, p r (n) is the center of the n th receiver and Z is the integer set. Denote λ as the Least Common Multiple (LCM) of {λ 1 , λ 2 , . . . , λ K }.
Proof. In a multi-path positioning model, Γ(k) in (16) is defined as: where Γ(k) is a block diagonal matrix: The n-th block Γ n (k) in the diagonal is a matrix with a size of M × L. The -th column of Γ n (k) is defined as e −iω k (τ n +τ ) , where τ n is an M × 1 column vector representing the propagation delays from the -th transponder to the M antennas in the n-th receiving station. Notice that: e −iω k (τ n +τ ) = e −iω k (τ n −τ n +τ n +τ ) = e −iω k (τ n −τ n ) e −iω k (τ n +τ ) a n (k)e −iω k (τ n +τ ) (A5) where τ n is the propagation delay from the -th transponder to the center of the n-th receiving array and a n (k) = e −iω k (τ n −τ n ) represents the array response of the n-th receiving station from the -th transponder.
If p e is a CSMC that satisfiesD 1,2 n = zλ and a 1 n (k) ≈ a 2 n (k), we get that e −iω k (τ 1 n +τ 1 ) ≈ e −iω k (τ 2 n +τ 2 ) from (A5). In this case, the 1 -th column of matrix Γ n (k) is approximately equal to the Notice that: w(k)α T =ᾱ T W(k), whereᾱ is defined in (38), and it is viewed as a reshaped form of the α. α is defined in (5), and: where I L is an identify matrix with a size of L × L and I N is an identify matrix with a size of N × N. Substitute (A9) into (A8): The eigenvalue decomposition of the second item of the right part is denoted by: where U =ŪΣ 1 2 . The right part of (A11) is simplified by: where: Sinceᾱ i+1 is the optimal solution of the relaxed programming, Q(η i+1 ) = H(ᾱ i+1 ) ≤ G(α i ,ᾱ i+1 , p e ) ≤ G(α i ,ᾱ i , p e ) = H(ᾱ i ) = Q(η i ),