Performance Analysis of the Direct Position Determination Method in the Presence of Array Model Errors

The direct position determination approach was recently presented as a promising technique for the localization of a transmitting source with accuracy higher than that of the conventional two-step localization method. In this paper, the theoretical performance of a direct position determination estimator proposed by Weiss is examined for situations in which the array model errors are present. Our study starts from a matrix eigen-perturbation result, which expresses the perturbation of eigenvalues as a function of the disturbance added to the Hermitian matrix. The first-order asymptotic expression of the positioning errors is presented, from which an analytical expression for the mean square error of the direct localization is available. Additionally, explicit formulas for computing the probabilities of a successful localization are deduced. Finally, Cramér–Rao bound expressions for the position estimation are derived for two cases: (1) array model errors are absent and (2) array model errors are present. The obtained Cramér-Rao bounds provide insights into the effects of the array model errors on the localization accuracy. Simulation results support and corroborate the theoretical developments made in this paper.


Introduction
The techniques of emitter localization using direction of arrival (DOA) measurements [1][2][3][4][5] play an important role in many areas, including vehicle navigation, localization and tracking of acoustic sources, and location services of satellite communications. In such localization systems, a single moving observer or multiple stationary observers are used to determinate the positions of the emitters. Generally, each observer is equipped with an antenna array for measuring the DOAs of the transmitted sources, and the emitter can then be located at the intersection of a set of lines of bearing [6][7][8]. The location procedure described above is typically called the two-step method. In the first step, the signal parameters (e.g., DOA [1][2][3][4][5], time difference of arrival (TDOA) [9,10], time of arrival (TOA) [11,12], frequency difference of arrival (FDOA) [13,14], frequency of arrival (FOA) [15], and received signal strength (RSS) [16,17]) are separately measured at several stations. In the second step, a central station uses the measurements to estimate the position coordinates of the sources. The two-step procedure is also known as the decentralized approach [18]. Note that although the two-step procedure is widely applied to the modern localization system, it is difficult to yield the optimal position estimate from the point of view of statistical characteristics. The reason is that the signal parameters are obtained by ignoring the constraint that all measurements must correspond to a common source position. As a result, information loss between the two steps is unavoidable. Although Table 1. Notational conventions.

Notation Explanation
⊗ Kronecker product Schur product diag [·] a diagonal matrix with diagonal entries formed from the vector blkdiag [·] a block-diagonal matrix formed from the matrices or vectors [·] † Moore-Penrose inverse of the matrix I n n × n identity matrix i (k) n the kth column vector of I n O n×m n × m matrix of zeros 1 n×1 n × 1 vector of ones λ max {·} the largest eigenvalue of the matrix || · || 2 Euclidean norm < · > n the nth entry of the vector < · > nm the nmth entry of the matrix Re{·} real part of the argument Im{·} imaginary part of the argument Pr{·} probability of the given event E [·] mathematical expectation of the random variable var [·] variance of the random variable

Time-Domain Signal Model
Consider an emitter and N base stations intercepting the transmitted signal. Each base station is equipped with an antenna array consisting of M elements. The transmitter's position is denoted by an L × 1 vector of coordinates p. In practice, L is equal to two or three, and cannot be larger than three. We consider the case where there is no multipath or non-line-of-sight (NLOS) phenomenon. The complex envelopes of the signal observed by the nth base station are then modeled by [29] x n (t) = β n a n (p)s(t − τ n (p ) − t 0 ) + ε n (t) (1 ≤ where • a n (p) is the nth array response to the signal transmitted from position p, • s(t − τ n (p ) − t 0 ) is the unknown signal waveform transmitted at unknown time t 0 , • τ n (p ) is the signal propagation time from the emitter to the nth base station (i.e., distance divided by signal propagation speed), • β n is an unknown complex scalar representing the channel attenuation between the transmitter and the nth base station, • ε n (t) is temporally white, circularly symmetric complex Gaussian random noise with zero mean and covariance matrix σ 2 ε I M . Assuming the observation vector x n (t) is sampled with period T, the kth sampled data can be expressed as x n,k = β n a n (p)s(kT − τ n (p ) − t 0 ) + ε n,k (1 ≤ k ≤ K) , where K is the number of snapshots.

Frequency-Domain Signal Model
To determinate the emitter position directly from all observations, it is desirable to separate the propagation delay τ n (p ) and transmit time t 0 from the signal waveform. This is easily achieved using the frequency-domain representation of the problem. Taking the discrete Fourier transform (DFT) of (2) produces [29] x n,k = β n a n (p)s k · exp{−jω k (τ n (p) + t 0 )} + ε n,k (1 ≤ n ≤ N ; 1 ≤ k ≤ K) , where • ω k = 2π(k − 1)/(KT) is the kth known discrete frequency point, • s k is the kth Fourier coefficient of the unknown signal corresponding to frequency ω k , • ε n,k is the kth Fourier coefficient of the random noise corresponding to frequency ω k .
It must be emphasized that the unknown and deterministic parameter set in (3) consists of p , t 0 , β n and s k . However, only the location vector p is of interest for the DPD approach. In addition, because the DFT is an orthogonal linear transformation, the distribution of the random noise vector ε n,k is the same as that of ε n,k , with first-and second-order moments given by. (1 ≤ k , l ≤ K ; k = l) .
Note that the DPD technique studied below is derived from (3).
It can be assumed, without loss of generality, that ||a n (p)|| 2 = ||s n || 2 = 1. Then, substituting (7) into (5) and applying algebraic manipulations leads to the concentrated problem [29] max p,s where A n (p) = blkdiag a n,1 (p) a n,2 (p) · · · a n,K (p) s = [s 1 · exp{−jω 1 t 0 } s 2 · exp{−jω 2 t 0 } · · · s K · exp{−jω K t 0 }] T (9) with a n,k (p) = a n (p) · exp{−jω k τ n (p)}. According to quadratic form theory, the cost function in (8) is maximized by selecting the vector s as the eigenvector corresponding to the largest eigenvalue of matrix N ∑ n=1 A H n (p)x n x H n A n (p). Therefore, (8) where B(p) = [A H 1 (p)x 1 A H 2 (p)x 2 · · · A H N (p)x N ] = A H (p)x (11) with A(p) = [A H 1 (p) A H 2 (p) · · · A H N (p)] H , It is important to stress that the second equality in (10) holds owing to the fact that given any matrix Z, the non-zero eigenvalues of Z H Z and ZZ H are identical [63]. Moreover, note that the dimensions of matrices B(p)B H (p) and B H (p)B(p) are respectively K × K and N × N. In practice, K is typically much greater than N and it is therefore more computationally efficient to perform the eigendecomposition on B H (p)B(p) instead of B(p)B H (p). Because the cost function in (10) is not a closed-form expression for p, the most straightforward method of solving (10) is to perform a grid search, as recommended in [29].
Note that when the location is estimated in multipath environments, the localization accuracy may obviously improve if the information contained in the non-line-of-sight signal components is exploited with the aid of appropriate channel modeling [35,36]. As a consequence, the signal model in (1) and (3) and the estimation criterion in (5) must be further adjusted to give a desired solution for the multipath model. Indeed, our performance analysis method also applies to the case of multipath propagation, but we only consider the single-path signal model in this paper owing to limited space.

Statistical Assumption and Effects of Array Model Errors
Assume that the actual array response, which differs from the nominal value, can be expressed aŝ a n (p) = a n (p) + ϕ n , (13) where ϕ n is the array model error. It must be emphasized that ϕ n is modeled as a stochastic variable rather than a deterministic variable throughout this paper. Moreover, there exist a variety of statistical assumptions that could be used to describe ϕ n in the literature. To make our results applicable to a more general situation, { ϕ n } 1≤n≤N is modeled as a set of independent complex Gaussian vectors with first-and second-order moments given by [43][44][45][46][47][48][49][50][51][52][53][54][55] Furthermore, array model error ϕ n is uncorrelated to sensor noise {ε n,k } 1≤k≤K for each base station. It is noteworthy that (14) will be used to determine the MSE and the CRB of the DPD estimator investigated in this paper.
When array model errors exist, the frequency-domain signal model in (3) becomeŝ where x n,k,0 is the true value ofx n,k in the absence of sensor noise and array model errors, and can be expressed as Defining the vectors and matrices it is easily verified from (15) and (17) that where From (17) and (18) we get where In the presence of array model errors, the emitter position is actually determined by whereB(p) = A H (p)X. We assume the optimal solution to (22) isp and its estimate error is p =p − p.
It is evident that the estimate error p depends on both sensor noise and array model errors. In where Obviously, ε c and ϕ c are related to sensor noise and array model errors, respectively. Further, we define two permutation matrices It can then be easily checked from (23) and (25) that ε c = Π ε ε * c and ϕ c = Π ϕ ϕ * c . In addition, it is straightforward to deduce from (17), (21), and (24) that E = O(||ε|| 2 ) and Ψ = O(|| ϕ|| 2 ).

MSE of Direct Position Determination Method in Presence of Array Model Errors
In this section, the MSE for the DPD method stated above is addressed in the presence of uncertainties in the model of the array manifold.

Perturbation Analysis on the Eigenvalues of Positive Semidefinite Matrix
Because the cost function in (22) is expressed as the maximal eigenvalue of some positive semidefinite matrix, an eigenvalue perturbation result is formally stated in a proposition as follows.

Proposition 1.
Let Z ∈ C N×N be a positive semidefinite matrix with eigenvalues λ 1 ≤ λ 2 ≤ · · · ≤ λ N , associated with unit eigenvectors u 1 , u 2 , · · · , u N , respectively. Moreover, λ n differs from the other eigenvalues. Assume Z is corrupted by a Hermitian error matrix Z ∈ C N×N , and the corresponding perturbed matrix is denotedẐ; i.e.,Ẑ = Z + Z ∈ C N×N . If the eigenvalues of matrixẐ are denotedλ 1 ≤λ 2 ≤ · · · ≤λ N , then the relationship betweenλ n and λ n can be described bŷ where The proof of Proposition 1 can be found in [21]. Note that Proposition 1 plays a fundamental role in our subsequent analysis.

Second-Order Perturbation Analysis on the Cost Function
Generally, first-order analysis is applied to predict the statistical performance of an estimator. The reason is that this analysis method gives the linear relationship between the estimation errors and measurement noise as well as model errors. As a result, the theoretical MSE of the estimator can be obtained according to statistical assumptions of the two sources of error. Moreover, first-order analysis is valid in most cases, provided that the error levels are not too high. In this paper, we employ this approach to derive the performance of the DPD estimator described above. For this purpose, first-order perturbation analysis is performed on the first derivative of the objective function in (22), or alternatively, second-order perturbation analysis is performed on the cost function in (22). Herein, because the analytical expression for the derivative of the cost function is rather complex, we prefer the second approach.
First, performing second-order perturbation analysis on matrixB(p) = A H (p)X leads tô where ξ = [ p T ε T ϕ T ] T consists of all error vectors, and The explicit expressions for .
A l (p) and ..
A l 1 l 2 (p) are given in Appendix A. It is seen from (28) and (29) that B (1) and B (2) collect all first-and second-order perturbation terms, respectively. It is deduced from (28) where From (31) and (32) we observe that C (1) and C (2) consist of all first-and second-order perturbation terms, respectively.
Let λ 1 ≤ λ 2 ≤ · · · ≤ λ N and u 1 , u 2 , · · · , u N be the eigenvalues and relevant unit eigenvectors of matrix C 0 , respectively. Additionally, it is not unreasonable to assume that the source location parameters are identifiable, which means that C 0 has unique maximal eigenvalue λ N . Meanwhile, it is noteworthy that the eigenvalue perturbation theory is extensively applied to the performance analysis in array signal processing for DOA estimation. To our best knowledge, there is no relevant mathematical tool that can be used to prove that the eigenvalues of C 0 are distinct. However, a large number of numerical investigations demonstrate that the possibility of the case of equal eigenvalues is small enough that we can ignore it. As a result, we define the matrix By combining Proposition 1 and (31), the cost-function value at pointp is given bŷ where It is seen from (35) that λ (1) N and λ (2) N group together all the first-and second-order error terms, respectively. The proof of (34) and (35) can be found in Appendix B. In the following, we express λ (1) N and λ (2) N as functions of ε c , ϕ c , and p. First, inserting the second equality in (29) into the first equality in (35) produces where in which { f k [· , ·]} 1≤k≤3 can be regarded as a set of vector functions, whose functional forms are given by The proof of (36) to (38) is provided in Appendix C. Secondly, substituting the second and third equalities in (29) into the second equality in (35) where in which {F ak [· , ·, ·]} 1≤k≤6 , {F bk [· , ·, ·]} 1≤k≤6 , {F ck [· , ·, ·]} 1≤k≤6 , and {G k [· , ·]} 1≤k≤3 can be viewed as matrix functions, which are given by The proof of (39) to (44) is provided in Appendix D. Substituting (36) and (39) back into (34) yieldŝ Evidently, Equation (45) can be considered as the second-order perturbation expression with respect to the error vectors ε c , ϕ c , and p. From (45), we get the linear relationship between the localization error p and sensor noise ε c as well as array model error ϕ c . The MSE of the DPD estimator can then be derived according to the statistical assumptions on the two sources of error.

MSE of Direct Position Determination Method
In light of the maximum principle, the true position p and estimated positionp satisfy the relations Obviously, the first equality in (46) leads to Additionally, using (45) and the second equality in (46), the localization error p =p − p is obtained by which further implies The second equality in (49) follows from (47). In (49), the linear relationship between the localization error p and the sensor noise ε c as well as the array model error ϕ c is formulated. It is easily observed from (49) that the positioning error vector p consists of two terms. The first term is associated with the sensor noise, which can be described as The second term is due to the array model errors, which can be written as According to the statistical assumptions in Sections 3 and 5, it is concluded that the localization error p is asymptotically Gaussian distributed with a zero mean and a covariance matrix given by where the second equality follows from (49) and the fact that ε c and ϕ c are statistically independent. Furthermore, (4), (14), and (23) together imply that Inserting (53) back into (52) leads to From (54) we see that the covariance matrix p is composed of two parts. The first part, due to the sensor noises, is expressed as The second part, due to the array model errors, is given by Remark 1. It is evident that the trace of P can be viewed as the MSE of localization errors under the combined effects of sensor noise and array model errors.

Remark 2.
When Φ (1) n → O and Φ (2) n → O , the trace of P can be viewed as the MSE of the localization errors when no array model errors are present. Moreover, the value of the trace of P approaches the CRB for the case of none of array model errors, which will be shown in Section 8.1. This is because the DPD method studied here is derived from the maximum likelihood (ML) criterion, which provides an asymptotically efficient solution.
Remark 3. When σ 2 ε → 0 , the trace of P can be used to quantify the sensitivity of positioning accuracy to array model errors, and represents the additional estimation errors resulting from uncertainties in the array manifold.

Remark 4.
It is easily seen from (55) and (56) that both P 1 and P 2 rely on matrix H 3 (p), which is the p-corner of the Hessian matrix of the cost function. If this matrix has a large condition number, the positioning accuracy might be high and, conversely, if this matrix is nearly singular, the location error may be extremely large. (54), it is observed that covariance matrix P is related to H 3 (p), H 5 (p), and H 6 (p). According to (38) and (40)-(44), the ijth element of matrix H 3 (p) is given by

Remark 5. From
In addition, the expressions for matrices H 5 (p) and H 6 (p) can be obtained from (38) and (40)- (44). However, the two formulas are complicated and we therefore have to omit them because of space limitations.

Success Probability of Direct Position Determination Method in Presence of Array Model Errors
The aim of this section is to deduce the success probability (SP) of the DPD method when array model errors exist. Two quantitative criterions are introduced to justify whether the localization is successful. Additionally, two analytical expressions for the SP of positioning are derived.

The First Success Probability of Direct Position Determination
It must be emphasized that the set of parameters {∆ l } 1≤l≤L in Definition 1 shall be appropriately chosen according to the practical scenario. The difference in these parameters reflects the importance of localization accuracy in distinct orientation. If the importance for each direction is identical, then these parameters can be set to the same value.
According to Definition 1, the joint probability density function of positioning error vector p is required for the calculation of the first localization SP. Applying the results in Section 6.3, the probability density function of random vector p is given by Consequently, the first localization SP can be determined by It is apparent from (59) that the first SP can be approximately obtained via numerical integration over a cube in high dimensional Euclidean space.
However, the high-dimensional numerical integration is not attractive from a computational viewpoint. If possible, it is preferable to get an explicit formula. Obviously, this is a non-trivial task and we only consider two-dimensional (2-D) localization scenarios (i.e., L = 2) for simplicity of mathematical analysis. First, an explicit formula with which to evaluate the joint probability of the Gaussian distribution is formally concluded in a proposition as below.
Proposition 2. Consider two joint Gaussian random variables z 1 and z 2 . The mean and variance of z 1 are m 1 and v 11 , respectively. The mean and variance of z 2 are m 2 and v 22 , respectively. In addition, the covariance of the two random variables is v 12 . It follows that where with Appendix E shows the proof of Proposition 2, which is along the lines of incomplete conditional moments theory presented in [46]. Note that Proposition 2 plays a significant role in the subsequent derivation process.
When L = 2, it can be verified by algebraic manipulation that The proof of (62) is shown in Appendix F. Applying the result in Proposition 2 and the definition of for arbitrary x is available from a table given in a textbook on probability theory.

Remark 7.
It must be pointed out that the above analytical results cannot be directly applied to the three-dimensional (3-D) case; i.e., L = 3. This can even be regarded as an open problem. Nevertheless, we can use numerical methods to compute this kind of SP in 3-D space. Indeed, there exist a number of efficient numerical integration methods with which to calculate the probability in (59), such as the Richardson extrapolation algorithm, Simpson algorithm, and Monte Carlo algorithm.

The Second Success Probability of Direct Position Determination
It is readily seen from Definition 2 that the second SP of positioning is equal to Pr {|| p|| 2 2 ≤ L∆ 2 }. To proceed, let us express p as p d = p 1/2 p 0 , where p 0 is a zero-mean Gaussian random vector with covariance matrix I L , and d = indicates that both sides have the same probability distribution. Consequently, || p|| 2 2 can be formulated as the quadratic form of p 0 : In light of the relationship between the cumulative distribution function and characteristic function [64], we have where φ || p|| 2 2 (t) denotes the characteristic function of || p|| 2 2 . Suppose that matrix P has eigenvalues γ 1 , γ 2 , · · · , γ L . Applying the property of the characteristic function, it can be proved that The substitution of (67) into (66) produces where Remark 8. It is clear from (68) that a one-dimension numerical integration over [0 , +∞) is required to evaluate the second SP. To this end, the values of the integrand shall be analyzed as t → 0 and t → +∞ .

Remark 9.
Applying L'Hospital's rule leads to Remark 10. The numerator of the integrand is bounded and the denominator tends to infinity when t → +∞ and, therefore, the integrand will be arbitrarily close to zero when t → +∞ . The integral upper limit in (68) can then be replaced by a sufficiently large positive number for the sake of simplicity.
Remark 11. It can be rigorously proved that the first SP is always smaller than the second SP, provided that ∆ 1 = · · · = ∆ L = ∆. The reason is that the first probability is computed by the numerical integral over a cube, while the second probability is equal to the integral over a circumscribed sphere of the cube.
As a byproduct of (68), we can present a new method of determining the radius of circular error probable (CEP), which is first defined in [65]. We denote r CEP by the radius of CEP, and it follows from its definition and (68) that As a consequence, a reasonable criterion for calculating r CEP is given by which can be solved via a one-dimensional grid search. In addition, it is noteworthy that although the solution for estimating r CEP is presented in [65], it is only applicable to 2-D localization scenarios. In contrast, the method proposed here is suitable for not only 2-D localization but also the 3-D scenario.

Cramér-Rao Bound on Covariance Matrix of Localization Errors
The CRB is a commonly used lower bound on the estimation error covariance of any unbiased estimator. In other words, the difference between the covariance and the CRB is a positive semi-definite matrix. Moreover, the CRB is expected to be a good predictor for the performance of the maximum likelihood estimator (MLE) at a moderate noise level. In this section, we derive the CRB for the estimate of the transmitter's position in two cases: (1) array model errors are absent and (2) array model errors are present. To this end, we first introduce the following proposition whose proof can be found in [66].

Proposition 3.
Assuming that the CRB matrix for the real vector z is equal to CRB(z), and defining a novel real vector as z = Jz, where J is an invertible matrix, the CRB matrix for vector z is given by CRB(z ) = J · CRB(z) · J T .

Cramér-Rao Bound on Position Estimate in Absence of Array Model Errors
This subsection is devoted to deriving the CRB for localization in the absence of array model errors. We begin by introducing a parameter vector that gathers all unknowns where To proceed, the data vector is defined as whose mean vector is given by Then, applying the results in [66,67], the CRB matrix for vector µ a can be obtained by where Using (16) and (77) and performing algebraic manipulations, the sub-matrices in (79) are described as where Note that only the p corner of the CRB matrix is of interest here. However, it is easily observed from (78) that matrix CRB(µ a ) does not exhibit a block-diagonal structure, because there might be correlation between the parameters. Hence, it is somewhat difficult to obtain the CRB for position vector p. To overcome this difficulty, we adopt the idea of [59,67] to redefine a parameter vector whose CRB matrix becomes block-diagonal. The new parameter vector is defined as where It is worth highlighting that because the vector µ a includes the source location parameters, it is meaningful to derive the CRB matrix for µ a . In addition, there is a one-to-one mapping between the new and old vectors µ and µ a . The relationship between them can be written in matrix form as where Then, combining the results in Proposition 3 and (84), the CRB matrix for µ a is given by where Combining (79), (83), and (87) leads to the orthogonal projection matrix where Inserting (88) back into (86) gives where We define three matrices The details of calculating the matrices in (92) are provided in Appendix G. Invoking the partitioned matrix inversion formula, the CRB matrix for position vector p is given by CRB(p) = σ 2 ε 2 · ((Re{V 1,1 }) −1 + (Re{V 1,1 }) −1 · Re{V 1,2 } · (Re{V 1,3 } − Re{V H 1,2 } · (Re{V 1,1 }) −1 · Re{V 1,2 }) −1 · Re{V H 1,2 } · (Re{V 1,1 }) −1 ) . (93)

Remark 12. The diagonal elements of CRB(p)
give the bounds for the estimation variance of the components in p when the array manifold is perfectly calibrated.

Remark 14.
Although there is no rigorous proof, it is expected that the trace of P 1 is asymptotically close to that of CRB(p). The reason for this is that the least square estimator in (5) is equivalent to the MLE, which is statistically efficient under the Gaussian noise model.

Remark 15.
By comparing the trace of CRB(p) with that of P, we can assess the expected degradation of the emitter location accuracy with respect to the amount of array model error. If the difference is significant, it can be concluded that the DPD method in [29] is sensitive to array model errors.

Cramér-Rao Bound on Position Estimate in Presence of Array Model Errors
This goal of this subsection is to derive the CRB for the position estimate in the presence of array uncertainties. Because in the present case the full parameter set contains both the deterministic parameters p, β, s, and σ 2 ε and the stochastic parameter ϕ, the CRB derivation should follow the Bayesian theory frame [68][69][70]. It is noteworthy that the CRB derivation can also be used for stochastic parameters, as processed in [68][69][70]. To this end, a novel parameter vector that comprises all the deterministic and stochastic unknowns is introduced where By performing similar algebraic manipulation in [68,69], the CRB matrix for vector µ b is formulated as where Note that (98) comes from the statistical assumption in (14). Appendix H provides the proof of (96).
Owing to the second term in the bracket of (96), it is impossible to get a CRB matrix with block diagonality as in (90) by linear transformation. As a result, the CRB matrix for position estimation can only be obtained from (96), although it may be computationally complex. Meanwhile, because the expressions for matrices Ω p , Ω Re{β} , Ω Im{β} , Ω Re{s} , and Ω Im{s} are given in (80), here we only need to deduce the expressions for matrices Ω Re{ ϕ} and Ω Im{ ϕ} . Applying (16) and (77) and performing algebraic manipulations gives The details of calculating the matrices in (101) appear in Appendix I. Through the application of the partitioned matrix inversion formula, the CRB matrix for position vector p is given by Note that the subscript "e" in (102) is used to distinguish the matrix CRB(p) for the case where the knowledge of the array manifold is accurate.

Remark 16. The trace of CRB e (p) is the bound for the localization MSE when array model errors exist.
Remark 17. It is apparent that the trace of CRB e (p) is larger than that of CRB(p) as the array model errors increase the uncertainties in parameter estimation.

Remark 18.
It can be readily proved that CRB(p) = CRB e (p) when Φ (1) n → O and Φ (2) n → O . Therefore, the CRB results derived in the presence of array model errors contain those for the case of no array model errors.

Remark 19.
Although there is no strict proof, it is not hard to conclude that the trace of P is greater than that of CRB e (p). The reason is that the DPD estimator discussed here does not take the array model errors into account and, thus, it is not statistically efficient for this case. Hence, a comparison of the trace of CRB e (p) with that of P allows us to decide whether a new DPD method that accounts for the array model errors is necessary to improve the emitter location accuracy.

Simulation Results
This section presents a set of Monte Carlo simulations to support the theoretical development in the previous sections. The empirical performances of the DPD method with and without array model errors are given, and they are compared both to the theoretical prediction values given in Sections 6 and 7 and to the CRBs presented in Section 8. The simulated values are averaged over 5000 independent trials. Moreover, the root-mean-square-error (RMSE), SP of localization, and radius of CEP are used to assess and compare the performance.

Discussion on RMSE of Direct Localization
This subsection focuses on the RMSE of the DPD method. Two sets of experiments are reported to illustrate the usefulness of the obtained results.

The First Set of Experiments
In the first set of experiments, the location estimation is performed on a 2-D plane and a simple array error model is used, which corresponds to case 1 in Section 4 in [44]. Specifically, ϕ n follows a circularly symmetric complex Gaussian distribution with second-order moments given by where σ ϕ is the standard deviation of the array model error.
The location geometry of the first set of experiments is shown in Figure 1

. The First Set of Experiments
In the first set of experiments, the location estimation is performed on a 2-D plane and a simple array error model is used, which corresponds to case 1 in Section IV in [44]. Specifically, n φ follows a circularly symmetric complex Gaussian distribution with second-order moments given by where φ  is the standard deviation of the array model error.
The location geometry of the first set of experiments is shown in Figure 1 (54) is in close agreement with the simulation result in the presence of array model errors. Consequently, the validity of the theoretical study in Section 6 is confirmed. Furthermore, when array model errors are absent, the empirical RMSE is very close to the CRB given by (92) and the theoretical RMSE in (55), which implies the asymptotical efficiency of the DPD method presented in [29], provided that the array is accurately calibrated. It is also seen that, as expected, the presence of array model errors leads to considerable deteriorations in location accuracy. Furthermore, Figures 2 and 6 show that the RMSE of the DPD method remains approximately constant no matter how much the SNR and sample number increase. The reason for this is that when the SNR or the sample number is large enough, the effects of sensor noise can be neglected and the localization errors are therefore primarily caused by array model errors, whose affects cannot be effectively eliminated in this DPD method yet. Additionally, we find that the RMSE performance in the presence of array uncertainties is significantly greater than the CRB provided by (102), especially when the standard deviation φ  increases (see Figure   3). Consequently, a new DPD method that accounts for array model errors is needed to improve the location accuracy.   (54) is in close agreement with the simulation result in the presence of array model errors. Consequently, the validity of the theoretical study in Section 6 is confirmed. Furthermore, when array model errors are absent, the empirical RMSE is very close to the CRB given by (92) and the theoretical RMSE in (55), which implies the asymptotical efficiency of the DPD method presented in [29], provided that the array is accurately calibrated. It is also seen that, as expected, the presence of array model errors leads to considerable deteriorations in location accuracy. Furthermore, Figures 2 and 6 show that the RMSE of the DPD method remains approximately constant no matter how much the SNR and sample number increase. The reason for this is that when the SNR or the sample number is large enough, the effects of sensor noise can be neglected and the localization errors are therefore primarily caused by array model errors, whose affects cannot be effectively eliminated in this DPD method yet. Additionally, we find that the RMSE performance in the presence of array uncertainties is significantly greater than the CRB provided by (102), especially when the standard deviation φ  increases (see Figure   3). Consequently, a new DPD method that accounts for array model errors is needed to improve the location accuracy.   (54) is in close agreement with the simulation result in the presence of array model errors. Consequently, the validity of the theoretical study in Section 6 is confirmed. Furthermore, when array model errors are absent, the empirical RMSE is very close to the CRB given by (92) and the theoretical RMSE in (55), which implies the asymptotical efficiency of the DPD method presented in [29], provided that the array is accurately calibrated. It is also seen that, as expected, the presence of array model errors leads to considerable deteriorations in location accuracy. Furthermore, Figures 2 and 6 show that the RMSE of the DPD method remains approximately constant no matter how much the SNR and sample number increase. The reason for this is that when the SNR or the sample number is large enough, the effects of sensor noise can be neglected and the localization errors are therefore primarily caused by array model errors, whose affects cannot be effectively eliminated in this DPD method yet. Additionally, we find that the RMSE performance in the presence of array uncertainties is significantly greater than the CRB provided by (102), especially when the standard deviation σ ϕ increases (see Figure 3). Consequently, a new DPD method that accounts for array model errors is needed to improve the location accuracy.

The Second Set of Experiments
In the second set of experiments, the source location is estimated in 3-D space and we assume that the array error is caused by sensor gain and phase uncertainties, which corresponds to case 2 in Section 4 in [44]. The second-order moments of ϕ n can then be expressed as where σ ϕ 1 and σ ϕ 2 are the sensor gain and phase perturbation standard deviation, respectively. Moreover, we assume σ ϕ 1 = 2σ ϕ 2 hereafter, and thus if σ ϕ 1 is changed, σ ϕ 2 alters accordingly.  Figures 8-12 show the RMSEs of the DPD method by varying the SNR of the emitter signal, standard deviation of sensor gain perturbation σ ϕ 1 , the number of array elements M, the ratio of array radius to wavelength, and the number of snapshots K.

The Second Set of Experiments
In the second set of experiments, the source location is estimated in 3-D space and we assume that the array error is caused by sensor gain and phase uncertainties, which corresponds to case 2 in Section IV in [44]. The second-order moments of n φ can then be expressed as  are the sensor gain and phase perturbation standard deviation, respectively. Moreover, we assume              Figure 11. RMSE of DPD as a function of ratio of array radius to wavelength. Figure 11. RMSE of DPD as a function of ratio of array radius to wavelength. Figure 11. RMSE of DPD as a function of ratio of array radius to wavelength. Owing to limited space, we do not present the results again. We simply highlight that the good agreement between the empirical and theoretical RMSE once again demonstrates the effectiveness of the theoretical development in Section 6.

Discussion on Success Probability of Direct Localization
This subsection focuses on the SP of the DPD method. Two sets of experiments are carried out to validate the obtained probability formulas, and the simulation parameters are the same as those in Section 9.1.

The First Set of Experiments
Both the localization scenario and the array error model for the first set of experiments are the same as those in Section 9.1.1. Moreover, the parameters 1 Δ and 2 Δ , which are used to specify the first SP, are set to the same value of 40, and the parameter Δ , which is related to the second SP, is also selected as 40. Because the localization scenario is on a 2-D plane, the theoretical value of the first SP can be obtained with (63). In Figures 13-15, we plot the two kinds of SP of the DPD method Owing to limited space, we do not present the results again. We simply highlight that the good agreement between the empirical and theoretical RMSE once again demonstrates the effectiveness of the theoretical development in Section 6.

Discussion on Success Probability of Direct Localization
This subsection focuses on the SP of the DPD method. Two sets of experiments are carried out to validate the obtained probability formulas, and the simulation parameters are the same as those in Section 9.1.

The First Set of Experiments
Both the localization scenario and the array error model for the first set of experiments are the same as those in Section 9.1.1. Moreover, the parameters ∆ 1 and ∆ 2 , which are used to specify the first SP, are set to the same value of 40, and the parameter ∆, which is related to the second SP, is also selected as 40. Because the localization scenario is on a 2-D plane, the theoretical value of the first SP can be obtained with (63). In Figures 13-15, we plot the two kinds of SP of the DPD method against the SNR of the emitter signal, standard deviation of the array model error σ ϕ , and number of snapshots K.    (63) and (68) can be supported. Furthermore, the simulated values in the absence of array model errors approach the lower bound calculated with the CRB in (93), which further indicates that the DPD estimator can achieve the CRB accuracy as long as the array is perfectly calibrated. However, when array model errors exist, the empirical values deviate significantly from the lower bound. Moreover, the difference increases with the standard deviation of array model error (see Figure 14). We thus need to develop a new DPD estimator with improved robustness against array model errors. Furthermore, it is seen that the first SP is always smaller than the second SP, which is consistent with the analysis in Remark 11.

The Second Set of Experiments
Both the localization scenario and the array error model for the second set of experiments are the same as those in Section 9.1.2. Because the situation is a 3-D localization scenario, the theoretical value of the first SP must be calculated with numerical integration methods. Herein, the Richardson extrapolation algorithm is exploited. Figures 16-18 (63) and (68) can be supported. Furthermore, the simulated values in the absence of array model errors approach the lower bound calculated with the CRB in (93), which further indicates that the DPD estimator can achieve the CRB accuracy as long as the array is perfectly calibrated. However, when array model errors exist, the empirical values deviate significantly from the lower bound. Moreover, the difference increases with the standard deviation of array model error (see Figure 14). We thus need to develop a new DPD estimator with improved robustness against array model errors. Furthermore, it is seen that the first SP is always smaller than the second SP, which is consistent with the analysis in Remark 11.

The Second Set of Experiments
Both the localization scenario and the array error model for the second set of experiments are the same as those in Section 9.1.2. Because the situation is a 3-D localization scenario, the theoretical value of the first SP must be calculated with numerical integration methods. Herein, the Richardson extrapolation algorithm is exploited. Figures 16-18 depict the two kinds of SP of the DPD method as functions of the SNR of the emitter signal, standard deviation of sensor gain perturbation σ ϕ 1 , and number of snapshots K. For Figures 16-18 we make similar observations as for Figures 13-15. We simply emphasize that the good agreement between empirical and analytical SP once again validates the probability formulas obtained in Section 7.

Discussion on Radius of CEP
This subsection discusses the radius of CEP of the DPD method. Two simulation experiments are conducted to illustrate the validity of (73), which is used to estimate the radius of CEP. The first and second simulation settings are the same as those in Figures 2 and 8, respectively. In the following two figures, the radius of CEP of the DPD method in the two experiments is plotted as a function of the SNR of the emitter signal. Figures 19 and 20 show that the simulation results agree well with the analytical results calculated with (73) and therefore the validity of (73) is corroborated. Moreover, we observe that the increase in the radius of CEP due to the array model errors is significant, especially when the SNR of the emitter signal is sufficiently high. Furthermore, when array model errors exist, the radius of CEP remains approximately constant no matter how much the SNR increases. Therefore, a robust DPD method that restrains the uncertainties in an array manifold is required. For Figures 16-18 we make similar observations as for Figures 13-15. We simply emphasize that the good agreement between empirical and analytical SP once again validates the probability formulas obtained in Section 7.

Discussion on Radius of CEP
This subsection discusses the radius of CEP of the DPD method. Two simulation experiments are conducted to illustrate the validity of (73), which is used to estimate the radius of CEP. The first and second simulation settings are the same as those in Figures 2 and 8, respectively. In the following two figures, the radius of CEP of the DPD method in the two experiments is plotted as a function of the SNR of the emitter signal.   Figures 19 and 20 show that the simulation results agree well with the analytical results calculated with (73) and therefore the validity of (73) is corroborated. Moreover, we observe that the  For Figures 16-18 we make similar observations as for Figures 13-15. We simply emphasize that the good agreement between empirical and analytical SP once again validates the probability formulas obtained in Section 7.

Discussion on Radius of CEP
This subsection discusses the radius of CEP of the DPD method. Two simulation experiments are conducted to illustrate the validity of (73), which is used to estimate the radius of CEP. The first and second simulation settings are the same as those in Figures 2 and 8, respectively. In the following two figures, the radius of CEP of the DPD method in the two experiments is plotted as a function of the SNR of the emitter signal.   Figures 19 and 20 show that the simulation results agree well with the analytical results calculated with (73) and therefore the validity of (73) is corroborated. Moreover, we observe that the

Conclusions
In this paper, the statistical performance of the DPD estimator presented in [29] is analytically studied when array model errors are present as well as signal waveforms are not available. The theoretical analysis begins with a matrix eigen-perturbation result, which can express the perturbation of eigenvalues as a function of the disturbance added to the Hermitian matrix. Then, the first-order asymptotic expression of the localization errors is given, from which the analytical expression for the MSE of the DPD estimator is obtained. Besides, the closed-form expressions for the calculation of the probabilities of a successful localization are also deduced, which can offer another theoretical perspective on the study of the localization accuracy. Additionally, the obtained probability formula can be used to provide a new criterion to estimate the radius of CEP. Finally, the CRB expressions for the position estimation are derived for two cases: (a) array model errors do not exist, and (b) array model errors are present and are drawn from Gaussian distribution. Several simulation experiments are performed to confirm the usefulness of the obtained results. The experimental results show that the uncertainties in the model of the array manifold can seriously deteriorate the source location accuracy of the DPD method. Therefore, our future work is to present a new DPD method that is expected to be more robust against the array model errors.

Appendix C. -Proof of (36) to (38)
From the second equality in (29), it follows for any vectors z 1 ∈ C K×1 and z 2 ∈ C N×1 that According to the last equality in (17), it can be readily checked that where f 1 [z 1 , z 2 ] is given in the first equality in (38). With (21), we have where f 2 [z 1 , z 2 ] is given in the second equality in (38). In addition, it can be easily verified that where f 3 [z 1 , z 2 ] is given in the third equality in (38). Combining (A9) to (A12) yields The substitution of (A13) into (A20) leads to where {F ck [· , ·, ·]} 1≤k≤6 are given in (43). For any vectors z 1 ∈ C K×1 and z 2 ∈ C N×1 , it can be easily verified from the third equality in (29) that According to the last equality in (17), we get where G 1 [z 1 , z 2 ] is given in the first equality in (44). It follows from (21) that where G 2 [z 1 , z 2 ] is given in the second equality in (44). In addition, it can be easily verified that < p > l 1 · < p > l 2 ·z H