Majorization-Minimization Method for Elliptic Localization in the Absence of Transmitter Position

This paper investigates the problem of elliptic localization in the absence of transmitter position. An efficient iterative method is developed to jointly evaluate the target and transmitter positions. Using the measurement information from the indirect paths reflected from the target and the direct paths between the transmitter and receivers, a non-convex maximum likelihood estimation (MLE) problem is formulated. Owing to the non-convex nature of the issue, we apply the majorization–minimization (MM) principle to address the MLE problem, which iteratively minimizes a convex surrogate function instead of the original objective function. Moreover, the proposed MM method is further extended to tackle a general scenario where both multiple unknown transmitters and receiver position errors are considered. Finally, numerical simulations demonstrate that the proposed MM method outperforms the state-of-the-art methods.


Introduction
Target localization has been a prevalent research topic with numerous applications in wireless sensor network (WSN) [1][2][3], communications [4], radars [5][6][7], and many others. A straightforward method to locate a target is direct localization [8], in which the position is determined directly from the received signals. The direct approach bears expensive computational complexity due to the multidimensional grid search. Hence, various types of indirect approaches have been presented based on different types of measurement information. The typical measurement information includes time of arrival (TOA) [9,10] , time-difference of arrival (TDOA) [11][12][13], received signal strength (RSS) [14], arrival of angle (AOA) [15], and their joint methods [16][17][18]. Among them, localization based on time information holds great research significance owing to its high-positioning accuracy. Elliptic localization is a time-based positioning model that is superior to some typical localization models, such as TOA-based and TDOA-based models [19]. In this model, the transmitter sends out a radar wave and several synchronous receivers collect the signal reflected from the interest target to determine its position. In addition, when the transmitter position is unknown, the signal received directly from the transmitter is also used to improve the positioning performance. Such a multistatic positioning approach [20] can offer more flexibility and better performance than the monostatic counterpart [21]. Hence, elliptic localization has been widely applied in multiple-input-multiple-output (MIMO) radars [22]. The flight time of the radar wave multiplied by its propagation speed produces the sum of the distances from the transmitter to the target and then to the receiver, which is referred to as the bistatic range [23,24]. The track of constant bistatic range in two-dimensional space is an ellipse with the transmitter and receiver positions as its foci, and the target position is on the ellipse. Thus, we can evaluate the target position by the intersection point of some ellipses produced by the bistatic ranges.
In reality, the transmitter position is potentially completely unavailable [20]. For instance, if the transmitter is set in some special inaccessible place or moves with time, its position cannot be accurately obtained. This also occurs in certain system designs such as the passive coherent location system [25][26][27], in which the transmitter position is purposely left unknown to simplify its structure and alleviate hardware requirements. Therefore, this situation bears important research significance. In [20], Zhang et al. first researched the elliptic localization problem in the absence of transmitter position. Using the measurements of the indirect and direct paths, a two-step weighted least squares (TSWLS) method was proposed for jointly estimating the target and transmitter positions. However, since the constrained relationships among those variables are difficult to directly and completely utilize, the TSWLS method has the threshold effect [28][29][30]. Zheng et al. [19] presented a semidefinite program (SDP) method to tackle the same problem, in which the semidefinite relaxation (SDR) technique is applied to convert this non-convex problem to a convex SDP problem. While the SDP method can approach the Cramér-Rao lower bound (CRLB) accuracy at a moderate noise level, this method requires a large number of complicated calculations arising from the CVX toolbox. Furthermore, there are several recursive methods that can be applied to solve the problem, such as the Gauss-Newton method [31] and the Quasi-Newton method [32,33].
This paper advances the topic by developing new solutions to the joint estimation problem. A maximum likelihood estimation (MLE) problem is formulated using the measurements of the indirect and direct paths. Since the problem is non-convex, its optimal solution cannot be solved directly. To tackle the difficulty, we propose an efficient iterative algorithm based on the majorization-minimization (MM) principle to address the formulated MLE problem. It requires constructing a convex surrogate function at each iteration that tightly upper bounds the original objective function, then minimizes the surrogate function to conduct the next iteration. Through the proposed algorithm, any given initial value will constantly approach the stationary point of the object function.
MM principle has been widely applied in signal processing, communications, and machine learning [34]. However, there is limited research about the MM principle for target localization. In fact, to the best of our knowledge, only Panwar et al. [35] has applied the MM principle to the multistatic target localization problem in terms of the precise known transmitter position. Furthermore, currently, no relevant works based on the MM principle have been discovered regarding the unknown transmitter position.
The main contributions of the current work are as follows: • We formulate an MLE problem that jointly estimates the target and transmitter positions, and propose an effective iterative algorithm based on the MM principle to solve the optimization problem. • We extend the proposed algorithm to a general scenario where both multiple transmitters at unknown positions and receiver position errors are considered. • We theoretically analyze the computational complexity and convergence of the proposed algorithm.
The rest of this paper is organized as follows. Section 2 provides the measurement model for elliptic localization. In Section 3, we present a brief overview of the MM technique followed by the derivation of the proposed approach. It is also utilized to solve the scenario where both multiple transmitters and receiver position errors are considered. We then conclude the section with a discussion on the complexity and convergence of the proposed algorithm. In Section 4, numerical simulations are presented and conclusions are given in Section 5.
Notations: Bold upper-case and bold lower-case letters denote matrices and vectors, respectively. The notations A T and A −1 represent the transpose and inverse of the matrix A, respectively. a i , a (i:j) , A ij and A (i:j,k:l) are the i-th element of a, the subvector containing the i-th to the j-th elements of a, the (i, j)-th element of A and the submatrix containing the elements from the i-th to the j-th row and the k-th to the l-th column of A, respectively. I p×p is the p × p identity matrix. O k×l is the k × l zero matrix. ⊗ and · represent the Kronecker product and the Euclidean norm, respectively.

Measurement Model
Let us consider an elliptic localization system for locating an unknown target, in which one transmitter and M receivers are deployed in p-dimensional space. As can be seen from Figure 1, the unknown transmitter represented by t o ∈ R p×1 sends out the radar signal. Several known receivers s j ∈ R p×1 for j = 1, · · · , M not only collect the indirect path signal reflected from the interest target denoted by u o ∈ R p×1 , but receive the direct path signal from the transmitter. and Multiplying by the signal propagation speed, the range measurements of the indirect and direct paths can be expressed by and respectively, where n r,j and n d,j are assumed to be independent zero-mean Gaussian variables with known variance σ 2 j and β 2 j , respectively.

Majorization-Minimization for Elliptic Localization
In this section, we present the proposed iterative algorithm based on the MM principle to tackle the elliptic location problem with an unknown transmitter position. Before we depict the proposed algorithm, we briefly introduce the MM principle.

Majorization-Minimization (MM)
The MM procedure is an efficient iterative method for addressing the non-convex problem. As shown in Figure 2, the MM procedure requires constructing a surrogate function (g(x|x k )) that tightly upper bounds the original objective function ( f (x)) at a given point (x k ), and then minimizes the surrogate function to obtain the next iteration. Figure 2. Illustration of MM principle.
As the tight upper bound of the f (x), the g(x|x k ) need satisfy Then in the minimization step, we update x as According to (3) and (4), we can simply obtain It can be readily proved from the inequality in (8) that the original objective function decreases monotonously during the updating process. The construction of the surrogate function is the key problem of the MM algorithm, which will determine the accuracy and complexity of the optimization approach. Interested readers can refer to [34] to acquire more details.

Elliptic Localization Using Single Transmitter
In this subsection, we consider the scenario with one transmitter and apply the MM framework to tackle the problem. From the measurement equations in (4) and (5), an MLE problem for evaluating the target position can be formulated as Combining u o and t o with θ, the cost function can be further rewritten as where θ = u oT , t oT T , Then, expanding the square of (10) and ignoring the constant terms, i.e., As can be observed, the problem is non-convex. Therefore, we can exploit the MM principle to develop an iterative method to solve it. As mentioned before, constructing a tight surrogate function is the key problem for the MM method [35]. Next, we shall derive the surrogate function by inequality scaling.
According to the inequality of arithmetic and geometric means, we known Plugging a = ( x x ) 2 and b = ( ỹ y ) 2 (x,x, y,ỹ are all positive numbers) into the above inequality, we can get Then we can easily obtain Using the above inequality, for any given θ k , the function ( D 1 θ + D 2 θ − s j ) 2 can be upper bounded as where According to the Cauchy-Schwartz inequality, for any given θ k , these functions D 1 θ , respectively, where Note that the denominators of these functions are not zero, and this situation can be avoided in practical scenario. Using these inequalities in (16) and (18)-(20), we get the following surrogate function for the cost function in (12) at (k + 1) th iteration: We can find that the surrogate function satisfies the conditions in (6), i.e., g(θ|θ k ) ≥ f (θ), ∀θ, g(θ k |θ k ) = f (θ k ). Therefore, it is one of the tight upper bounds of the original cost function ( f (θ)) at the point (θ k ). Then, minimizing the surrogate function, the updating process is expressed as It can be simply seen that the surrogate function is convex. By taking the derivative of (22) and letting it equal zero, we can get where A limit on the number of iterations, K max , should be applied for practical scenarios. The criterion to stop the iteration is when or when the time of iteration attains K max , where is a preset threshold. Then, the estimated results of the target and transmitter positions are expressed as respectively, where θ * is the final result of iterations. In Algorithm 1, we give the relevant pseudo code.

Elliptic Localization Using Multiple Transmitters with Receiver Position Errors
It is common in practical applications, especially in sonar/radar, to deploy multiple transmitters for improving performance. Moreover, due to imperfections in GPS precision, the available receiver positions may present errors. Ignoring these errors can generate significant performance degradation [20]. In the subsection, we further apply the MM method to tackle the general scenario where both multiple unknown transmitters and receiver position errors are considered.
The unknown position of the i-th transmitter is represented by t o i ∈ R p×1 , i = 1, · · · , N. Then, the range measurements of the indirect and direct paths are expressed by and respectively. Furthermore, lets j ∈ R p×1 , j = 1, · · · , M be the available position of the j-th receiver. The true position s o j is not known and always modeled as where ∆s j is the position error of s j that obeys the zero-mean Gaussian distribution with known covariance γ 2 j I p×p . From the measurement equations in (28)- (30), an MLE problem for evaluating the target position can be formulated as Using a vector θ to combine these unknown variables, we can obtain where, Then, expanding the square of (32) and ignoring the constant terms, i.e., Similarly, we shall utilize the MM principle to account for the non-convex problem as well. As already mentioned, we shall utilize these inequalities mentioned in the previous subsection to construct the surrogate function. Using the inequality in (15), for any given θ k , the function ( H i θ + P j θ ) 2 can be upper bounded as According to the Cauchy-Schwartz inequality, for any given θ k , these functions H i θ , P j θ and Q ij θ can be upper bounded as and respectively, wherex Note that the denominators of these functions are not zero as well.
Using these inequalities in (35) and (37)-(39), we get the following surrogate function for the cost function in (34) at (k + 1) th iteration: Minimizing the surrogate function, the updating process can be expressed as where The criterion to stop the iteration is same to the previous subsection. Then, the estimated results of the target, transmitter, and receiver positions are expressed as respectively, where θ * is the final result of iterations. As we can see from (44), different from the existing approaches in [19,20], the proposed method can not only estimate the positions of target and transmitters, but also calibrate the erroneous receiver positions.

Computational Complexity and Convergence of the Proposed Algorithm
We first analyze the computational complexity of the proposed methods. In each iteration, the complexity is mainly dominated by computing (24) for the single transmitter case and (42) for the multiple transmitters case. Hence, their complexities are on the order of O(M 2 ) and O((N + M) 4 ), respectively. Moreover, the average computation time will be shown in the next section.
Next, we shall discuss the convergence of the proposed algorithm. Since the algorithm is based on the MM framework, it is seen from (8) that the objective function monotonically decreases with the iteration. Moreover, the surrogate function at each iteration is convex, so we can easily get its solution θ k+1 in (k + 1) th iteration that makes and g(θ k+1 |θ k ) ≤ g(θ|θ k ).
Hence, through the proposed algorithm, any given initial value will constantly approach the stationary point of the object function.

Simulation Results
In this section, numerical simulations are conducted to evaluate the performance of the proposed method in 2-D for ease of illustration. The proposed method is compared with the SDP method [19] (labeled as "SDP"), the two-stage WLS method [20] (labeled as "TSWLS") and the Gauss-Newton method (labeled as "Gauss-Newton", the derivation is shown in Appendix A). Moreover, the CRLB is also included as a benchmark. Owing to many unknown variables in the objective function, we choose to initialize θ by successive estimation. The successive estimation can be regarded as a TOA localization problem (based on direct paths) and an elliptic localization problem (based on indirect paths). The TOA localization problem is solved to evaluate the transmitter position by the LS method in [36]. Then, using the evaluated transmitter position, the elliptic localization problem is solved to evaluate the target position by the LS method in [28]. Meanwhile, the receiver position errors are not considered. The noise variances of the indirect and direct paths are set to the same value (represented by σ 2 ). The maximum number of iterations is set to K max = 4000 and the predetermined threshold is set as = 10 −3 .
The root mean square error (RMSE) is utilized to evaluate the accuracy of these algorithms, which is expressed as where G are the number of geometry configurations and L is the times of Monte Carlo (MC) runs for each geometry configuration. η * gl is the estimated result in the l-th MC run for the g-th configuration. G and L are adopted as 10 and 500 in the simulations.
Next, we shall test those algorithms in the single transmitter and multiple transmitters cases, respectively.  Figure 3. We find that the RMSE of the TSWLS method attains the CRLB accuracy when 10log 10 (σ 2 ) ≤ 5 dBm 2 , whereas deviates from the CRLB when 10log 10 (σ 2 ) ≥ 10 dBm 2 . This may be due to the threshold effect of the TSWLS method. By contrast, the proposed method, SDP method and Gauss-Newton method show better threshold behavior, and the first has a higher performance because its algorithm is based on a more accurate problem formulation. Figure 4 displays the estimated results of the transmitter position. It is observed that the RMSEs of the proposed method, SDP method and Gauss-Newton method attain the optimal performance when 10log 10 (σ 2 ) ≤ 35 dBm 2 , and slightly deviates from the CRLB at 10log 10 (σ 2 ) = 40 dBm 2 . Scenario 2: In this scenario, receiver position errors are considered to be present. Assuming the covariance matrix of the receiver position errors Q s = σ 2 s J, where σ 2 s relates to the sensor position error power and J = diag [5,20,15,10] ⊗ I 2×2 . The deployment of the target, transmitter, and receiver is the same as that in Scenario 1. Keeping σ 2 s at 1 m 2 , Figure 5 confirms that the RMSE of the TSWLS method attains the CRLB accuracy when 10log 10 (σ 2 ) ≤ 5 dBm 2 and deviates from the CRLB when 10log 10 (σ 2 ) ≥ 10 dBm 2 much earlier than the other methods. When 10log 10 (σ 2 ) = 40 dBm 2 , the proposed methods still performs better than the SDP method and Gauss-Newton method.   [5,20,35,10,25,30,15,40] ⊗ I 2×2 . The deployment of the actual coordinates of the target, transmitter, and receiver is the same as that in Scenario 1. As can be seen in Figure 6, the performance of these methods improves as the number of receivers increases. The RMSE of the TSWLS method deviates from the CRLB accuracy when M ≤ 4. The proposed method, SDP method and Gauss-Newton method attain the CRLB accuracy at the entire tested range except for M = 3. For M = 3, the SDP method obviously deviates from CRLB; however, the proposed method and Gauss-Newton method just slightly deviate from CRLB, and the former has better performance.

Multiple Transmitters Case
Scenario 4: Multiple transmitters and receiver position errors are simultaneously considered in this Scenario. Three transmitters and four receivers are used, and the other parameter settings are the same as that in Scenario 2. As can be seen in Figure 7, the observations are similar to those in Figure 5, whereas the RMSEs are lower at the same σ 2 . This is caused by using multiple transmitters. Scenario 5: Fixing three transmitters, the number of receiver M increases from 3 to 8. We let σ 2 = 100 m 2 and Q s = σ 2 s J M with σ 2 s = 1 m 2 , and the other parameter settings are the same as that in Scenario 3. As we can see in Figure 8, the performance of these methods improves with the rise of the receiver number. Moreover, the TSWLS method deviates from the CRLB accuracy when M ≤ 4 and the SDP method and Gauss-Newton method deviate from the CRLB accuracy at M = 3. By contrast, the proposed method can approach the CRLB accuracy at the whole tested noise level. Finally, the average CPU time of the methods in Scenario 1 and Scenario 4 is shown in Table 1. The simulation data is obtained by utilizing a PC with Intel Core i7 3.2 GHz processor. As reported in Table 1, the computational complexity of these methods increases as the measurement model becomes more complex. Furthermore, the computational complexity of the proposed method is higher than that of the TSWLS method and Gauss-Newton method, and lower than that of the SDP method.

Conclusions
In this paper, we investigate an elliptic localization problem when the transmitter position is unavailable. An MM method is developed to jointly estimate the target and transmitter positions. Subsequently, the MM method is further extended to account for scenarios where both multiple unknown transmitters and receiver position errors are considered. Simulation results demonstrate that the proposed method outperforms stateof-the-art methods and the erroneous receiver positions can also be calibrated.

Abbreviations
The following abbreviations are used in this manuscript: Finally, the updating of θ is Then (A10) can be recast as Finally, the updating of θ is