Real-Valued Direct Position Determination of Quasi-Stationary Signals for Nested Arrays: Khatri–Rao Subspace and Unitary Transformation

The features of quasi-stationary signals (QSS) are considered to be in a direct position determination (DPD) framework, and a real-valued DPD algorithm of QSS for nested arrays is proposed. By stacking the vectorization form of the signal’s covariance for different frames and further eliminating noise, a new noise-eliminated received signal matrix is obtained first. Then, the combination of the Khatri–Rao subspace method and subspace data fusion method was performed to form the cost function. High complexity can be reduced by matrix reconstruction, including the modification of the dimension-reduced matrix and unitary transformation. Ultimately, the advantage of lower complexity, compared with the previous algorithm, is verified by complexity analysis, and the superiority over the existing algorithms, in terms of the maximum number of identifiable sources, estimation accuracy, and resolution, are corroborated by some simulation results.


Introduction
Source localization technology is an essential part of many fields, including rescue operation, resource exploration, intelligent transportation, and underwater detection [1][2][3][4]. Initially, typical localization methods, such as time of arrival (TOA) [5,6], angle of arrival (AOA) [7,8], and frequency difference of arrival (FDOA) [9,10] are always performed in a two-step mechanism. The intermediate parameters containing information regarding the source position are estimated first, such that the source position can be determined by methods that are based on the geometric relationship between the parameters. Nevertheless, the two-step algorithm cannot perform optimally, due to the inevitable information loss between the two steps. Moreover, a sharp decline in accuracy will be caused if parameter matching errors occur in multiple source scenarios.
To circumvent these problems in the two-step framework, a new one, called direct position determination (DPD), is first proposed in [11]. As its name suggests, the DPD algorithm is performed by processing the raw received signals to determine the source position. Thus, it skips the step of intermediate parameter estimation and takes the correlation among different received signals into account. The research results in [11] show that much higher accuracy can be achieved by DPD algorithms, compared with two-step algorithms, especially under low signal-to-noise ratio (SNR) conditions.
As the information associated with the source position cannot be completely ignored in DPD, a series of DPD algorithms based on different information types have been proposed. The information regarding TOA and AOA was considered in [11], where the maximum 1.
The features of QSS are considered in the DPD model for NAs, where the cost function is constructed by combining the Khatri-Rao subspace and SDF methods, and the QSS-SDF-DPD algorithm is derived.

2.
The original dimension-reduced matrix is modified, and the unitary transformation method is exploited for the purpose of releasing the computational burden; then, the R-QSS-SDF-DPD algorithm is proposed. 3.
Apart from the given Cramer-Rao bound (CRB), the complexity analysis, summary of advantages, superiority of the R-QSS-SDF-DPD algorithm in terms of computational complexity, maximum number of identifiable sources, localization accuracy, and sources resolution is confirmed by some simulation experiments.
The remaining parts of this paper are organized as follows. Section 2 presents the system model of QSS-DPD for NAs. The proposed algorithm, which contains the Khatri-Rao subspace method for NA and matrix reconstruction for complexity reduction, is derived in Section 3, and a summary is given at the end. Section 4 provides the CRB, complexity analysis, and advantages of the proposed algorithm. Then, some relevant simulation results are shown in Section 5. The last section draws some conclusions on the paper. Notation: Throughout the paper, the upper-case bold character and upper-case one are used to represent a matrix and vector, respectively, and a variable is denoted with a lower-case character. |·| and · represent the operation of taking the magnitude and Euclidean norm of a vector, respectively. E{·}, vec(·), diag(·), blkdiag{·}, and Re(·) represent the operation of expectation, vectorization, diagonalization, block diagonalization, and taking the real part, respectively. ⊗ and represent the operation of the Kronecker and Khatri-Rao products, respectively. C M×N and R M×N represent the complex number set and real number set with dimension M × N, respectively. The operation of inverse, conjugate, transpose, and conjugate transpose are represented by (·) −1 , (·) * , (·) T , and (·) H , respectively. arcsin(·) represents the operation of arcsine. I n , J n , 0 N , and 1 n represent the n × n identity matrix, row-flipped form of the n × n identity matrix, n × 1 vector with all zeros, and n × 1 vector with all ones.

System Model
Consider the two-dimension scenario presented in Figure 1, where the K (it is assumed to be known, as it can be estimated by some methods [32][33][34][35]) far-field narrowband uncorrelated sources are intercepted by N base stations, which are equipped with a NA. As the location of the base stations are known, assume they are located at v n = [x v n , y v n ] (n = 1, 2, . . . , N), and the sources are located at q k = [x k , y k ] (k = 1, 2, . . . , K). The specific structure of the M element NA that is exploited in this scenario is shown in Figure 2, where the first and second levels consist of M 1 and M 2 elements (M = M 1 + M 2 ), respectively. The place of all physical array elements can be included in a set Θ, given by [21]: where d denotes the unit adjacent spacing.  Assume the kth source is impinging on the nth base station from θ n,k ; then, the received signal vector intercepted by the nth base station at the tth (t = 1, 2, . . . , T) sampling time can be presented by [14]: x n (t) = K ∑ k=1 a n (θ n,k )s n,k (t) + n n (t) (2) where a n (θ n,k ) = e −j2πd 1 sin θ n,k /λ , . . . , e −j2πd m sin θ n,k /λ , . . . , e −j2πd M sin θ n,k /λ T denotes the steering vector of θ n,k = arcsin((x k − x v n )/ q k − v n ), which is the DOA from the kth source to the nth base station, d m ∈ Θ (m = 1, 2, . . . , M), λ is the signal wavelength, s n,k (t) represents the envelope of kth source incident on the nth base station, and n l (t) is the Gaussian white noise vector.
According to Equations (8) and (9), a local covariance matrix can be defined by [27]: where Λ n, f = diag g n,1, f , g n,2, f , . . . , g n,K, f denotes the local source covariance matrix of the f th frame, and C n is the spatial noise covariance. According to Property 1 in [27], the vectorization of R n, f in Equation (10) can be expressed as: where g n, f = g n,1, f , g n,2, f , . . . , g n,K, f T ∈ R K×1 is a real-value vector consisting of the diagonal elements in Λ n, f .
Interestingly, a new received signal vector y n, f is obtained, whose steering matrix is A * n A n and noise vector is vec(C n ). By stacking the vectorization forms of all local covariance matrices, the new received signal matrix Y n ∈ C M 2 ×F is obtained: where G n = g n,1 , g n,2 , . . . , g n,F ∈ R K×F is the new signal matrix.
Based on Assumption A4 in [27], the covariance of the source signal is time-varying. Therefore, G n can be treated as a new incoherent source signal matrix with F snapshots. Note that this is different from the general case, where the virtualization of the received signal makes the signals coherent, and some decoherence operations should be performed. In contrast, that problem does not occur in this paper, and the advantages of high DOF and large aperture can be fully preserved. As shown in Figure 3, we compare the DOF of the general signal model (the spatial smoothing method in [21] is adopted to overcome the coherent signal problem that is caused by the virtualization) and QSS model in this paper. Obviously, the QSS model has approximately twice as many DOF as the general signal model.
when postmultiplying Y n by it, the unknown noise in Equation (12) can be eliminated, and the new noise-eliminated received signal matrixỸ n ∈ C M 2 ×F can be given by [27]: According to Assumption A5 in [27], G n is a full row-rank matrix if F ≥ K + 1, which means rank G T n = K. Note that H ⊥ is non-singular, and G n H ⊥ is a full row-rank matrix. The noise subspace can be obtained after performing singular value decomposition (SVD) onỸ n or eigenvalue decomposition (EVD) on its covariance matrix, so that the SDF algorithm [14] can be exploited to estimate the source positions.
However, some direct processing onỸ n is time-consuming. To avoid this problem, the R-QSS-SDF-DPD will be proposed in the next section.

The Proposed Algorithm
In this section, the derivation process of the QSS-SDF-DPD and its modified realvalued version (R-QSS-SDF-DPD) are given in detail, including the Khatri-Rao subspace method for NA and matrix reconstruction for complexity reduction.

Khatri-Rao Subspace Method for NA
For a uniform linear array (ULA) with L elements, the new steering matrix A * n A n can be characterized as [27]: where: ×K denotes the dimension-reduced virtual array steering matrix, and b(θ k ) = [e j2π(L−1)d sin θ k /λ , . . . , e j2πd sin θ k /λ , 1, e −j2πd sin θ k /λ , . . . , e −j2π(L−1)d sin θ k /λ ] T denotes the virtual steering vector of θ k , which is the DOA of the kth source. However, due to the discontinuity between the two-level elements of the NA, Equation (14) cannot be directly applied to it. To overcome this difficulty, a transformation should be applied to the original array steering matrix A n .
For the two-level NA exploited in this paper, a total of 2R N A − 1 (R N A = M 2 (M 1 + 1)) consecutive distinct elements can be obtained after the operation of Khatri-Rao product A * n A n . If we take the non-negative part of them as the virtual ULA, then the original NA can be treated as a subset of the virtual ULA, and their array steering matrices scarify [29]: where vector with 1 at the jth (j = 1 + d m /d) entry and 0 elsewhere [29].B n = b n,1 ,b n,2 , . . . ,b n,K ∈ C R N A ×K represents the array steering matrix of the virtual ULA, andb n,k = [1, e −j2πd sin θ n,k /λ , . . . , e −j2π(R N A −1)d sin θ n,k /λ ] T ∈ C R N A ×1 represents the corresponding steering vector of the kth source. Substitute Equations (16) and (14) into Equation (13); then,Ỹ n can be rewritten as [29]: is the dimension-reduced matrix with the form of Equation (15),B n = b n (θ n,1 ),b n (θ n,2 ), . . . ,b n (θ n,K ) ∈ C (2R N A −1)×K denotes the dimension-reduced virtual array steering matrix with the form of B in Equation (14), and its kth column vector b n (θ n,k ) ∈ C (2R N A −1)×1 represents the corresponding steering vector with the form of b(θ k ) in B.
According to the definition ofΓ and P, they are both column orthogonal, so (P ⊗ P)Γ is also column orthogonal, which is easy to verify [29].
, and the dimension ofỸ n can be reduced after a linear transformation, which is given by [29]: where Y n is the dimension-reduced, noise-eliminated received signal matrix, andW = V −1/2 (P ⊗ P)Γ T is the dimension-reduced matrix. As presented in Equation (18), the dimension of the original received signal matrix is reduced from R 2 N A to 2R N A − 1. It can be seen from Equation (18) that a virtual noise-eliminated received signal matrix Y n is obtained, whose array steering matrix isṼ 1/2B n and source signal matrix is G n H ⊥ ∈ C K×F . Therefore, the covariance matrix of the dimension-reduced, noiseeliminated received signal matrix R n ∈ C (2R N A −1)×(2R N A −1) can be calculated by: After the eigenvalue decomposition ofR n , we obtain the corresponding signal sub-spaceÊ S n and noise subspaceÊ N n , which are given by: whereΛ S n = diag σ 2 n,1 , σ 2 n,2 , . . . , σ 2 n,K andΛ N n = diag σ 2 n,K+1 , σ 2 n,K+2 , . . . , σ 2 n,2R N A −1 σ 2 n,1 > σ 2 n,2 , . . . , σ 2 n,K > σ 2 n,K+1 ≥ . . . ≥ σ 2 n,2R N A −1 denote the diagonal matrices made up of the maximum K eigenvalues and remaining ones, respectively. The corresponding eigenvectors formÊ S n andÊ N n , respectively. Hence, according to the SDF algorithm in [14], the spectrum function of the QSS-SDF-DPD algorithm for NAs sp QSS−SDF (p) can be constructed by: where p is the position of search grid point in the search area Ξ, andb n (θ(p)) ∈ C (2R N A −1)×1 represents the steering vector of p for the nth base station, which has the same form as b(θ k ) in B. After finding the maximum K points of Equation (21), all source positions can be estimated.

Matrix Reconstruction for Complexity Reduction
Note that high complexity cannot be avoided during the peak search process of sp SDF (p), due to the existence ofṼ 1/2 and lots of complex-valued multiplications. For the purpose of releasing the computational burden, the dimension-reduced matrix will be modified, and the unitary transformation method [30,31] will be adopted.
It can be observed from Equation (18) that, if we pre-multiply Y n byṼ −1/2 , we can obtain: where Y M n is the modified, dimension-reduced, noise-eliminated received signal matrix, andW M =Ṽ −1 (P ⊗ P)Γ T ∈ R (2R N A −1)×M 2 is the modified dimension-reduced matrix.
Interestingly and coincidentally, for the NA,W M is equal to the Equation (9) in [39], which is expressed by: where ∆ j,k = d j − d k /d d j , d k ∈ Θ represents the position difference between d j and d k , and ω ∆ j,k represents the number of pairs d j , d k , whose difference d j − d k /d is equal to ∆ j,k .
According to Equation (22), Y M n can be treated as a modified noise-eliminated received signal matrix, whose array steering matrix isB n and source signal matrix is G n H ⊥ ∈ R K×F . Note that, compared with the virtual array steering matrixṼ 1/2B n before modification, V 1/2 is eliminated, which means the matrixṼ 1/2 would be removed from sp SDF (p). Thus, it partially reduces the computational burden.
As the virtual array configuration formingB n is centrosymmetric, the unitary transformation method [30,31] can be directly adopted.
Define the unitary matrix by [30]: Then, a real-valued matrix Y MR n ∈ R (2R N A −1)×F can be obtained by the unitary transformation given by: Correspondingly, a real-valued covariance matrixR Hence, after performing eigenvalue decomposition onR Finally, fuse all the real-value noise subspace of all base stations, and the spectrum function of R-QSS-SDF-DPD algorithm for NAs sp R−QSS−SDF (p) can be expressed by: is the real-valued peak search steering vector. Then, find the maximum K points of sp MR−QSS−SDF (p), and take the corresponding coordinate positions as the estimates of all source positions.

Summary of The Proposed Algorithm
As a summary, the main steps of the QSS-SDF-DPD and R-QSS-SDF-DPD algorithms are listed below, where the first six steps belong to the former, and the last four ones belong to the latter.
Step 1: Calculate the covariance matrix of the received signal for each frame R n, f by Equation (10).
Step 2: Construct the vectorization forms of all R n, f , and stack them together to obtain the new received signal matrix Y n by Equation (12).
Step 4: Obtain the original dimension-reduced, noise-eliminated received signal matrix Y n by Equation (18).
Step 5: Calculate the covariance matrix of Y n by Equation (19), and obtain the noise subspace matrixÊ N n by Equation (20).
Step 6: Construct the original spectrum function of QSS-SDF-DPD sp QSS−SDF (p) by Equation (21), and estimate all source positions by finding the maximum K points.
Step 7: Modify the dimension-reduced matrix to obtain the modified noise-eliminated received signal matrix Y M n according to Equation (22).
Step 8: Obtain the real-valued matrix Y MR n , according to Equation (25).
Step 9: Calculate the real-valued covariance matrixR MR n by Equation (26), and obtain the real-valued noise subspace matrixÊ MRN n by Equation (27).
Step 10: Construct the spectrum function of R-QSS-SDF-DPD sp R−QSS−SDF (p) by Equation (28), and estimate all source positions by finding the maximum K points.

CRB
Based on the derivation results in [24,40], the CRB of DPD for NAs can be expressed by: where σ 2 is the variance of the noise, and: denote the partial derivatives of A, with respect to x k and y k , respectively.

Complexity Analysis
In this subsection, the complexity of the SDF-DPD, QSS-SDF-DPD, and R-QSS-SDF-DPD algorithms are compared. For the sake of fairness, the two-level NA is employed in all algorithms, and the spatial smoothing method [21] is adopted in SDF-DPD algorithm. Moreover, we count the number of real-valued multiplication operations, instead of the complex-valued multiplication operations, which are equivalent to four real-valued multiplication operations. Table 1 presents the results of the comparison, and the corresponding intuitive form is depicted in Figure 4, where the search area is assumed to be divided into N x × N y grids, M = M 1 + M 2 , and R N A = M 2 (M 1 + 1).
It can be concluded that the proposed R-QSS-SDF-DPD algorithm has a much lower complexity than the QSS-SDF-DPD algorithm, and it is slightly higher than the SDF-DPD algorithm, which confirms the effectiveness of the dimension-reduced matrix modification and unitary transformation. This means that the R-QSS-SDF-DPD algorithm is more practical in engineering applications.

Algorithm Major Steps Complexity (Real-Valued Multiplication Operation)
The SDF-DPD The proposed QSS-SDF-DPD The proposed R-QSS-SDF-DPD

Advantages
Due to the utilization of the QSS features, modification of the dimension-reduced matrix, and unitary transformation, the proposed R-QSS-SDF-DPD algorithm possesses the following advantages, when compared to the existing algorithms.

1.
More sources can be estimated than the traditional SDF-DPD algorithm, even when K > R N A − 1 (K is the actual number of sources, and R N A − 1 is the maximum number of identifiable sources for two-level NA exploiting spatial smoothing method); 2.
Larger array aperture, lower localization error, and higher resolution can be obtained, compared to the SDF-DPD algorithm; 3.
Less computational burden than the QSS-SDF-DPD algorithm, before dimensionreduced matrix modification and unitary transformation.

Simulation Results
We compare the estimated performance of the SDF-DPD, QSS-SDF-DPD, and R-QSS-SDF-DPD algorithms by calculating the root mean square error (RMSE), as defined by: where N E denotes the number of Monte Carlo simulation experiments, and [x k,ne ,ŷ k,ne ] denotes the estimate of [x k , y k ] in the neth simulation experiment. N E is set to be 1000 in all the following simulation experiments, unless a special statement is given. Besides, in order to compare the resolution of algorithms, the resolving probability P r [41] is defined by: where N s is the number of simulation experiments in which the two sources are successfully distinguished. In this section, two sources are placed parallel to the x-axis, one of them is fixed at q 1 = [x 1 , y 1 ]m, while the abscissa of the other q 2 = [x 1 + ∆x, y 2 ]m is adjusted, where ∆x denotes the distance of the two sources. For each simulation experiment, if the locations of two sources are estimated and satisfy |x 1 − x 1 | < ∆x/2, |x 2 − x 2 | < ∆x/2 (x 1 andx 2 are the estimated abscissa of the two sources), then this simulation experiment passes; otherwise, it fails. Simulation 1: To verify the feasibility of the proposed R-QSS-SDF-DPD algorithm, the scatter plot is shown in Figure 5. In this simulation, the number of sources is K = 7, and they are located at q 1 Figure 5, the proposed algorithm can estimate all source positions, even when K = 7 > R N A − 1 = 5, which proves the first advantage that was mentioned in Section 4.3. Besides, the feasibility of the R-QSS-SDF-DPD algorithm, in regard to estimate accuracy, is also confirmed. The structure of array is M = 4 for ULA, and M 1 = 2, M 2 = 2 for NA. The number of frames is F = 80, each frame consists of T F = 500 snapshots, and the SNR varies from −6 dB to 18 dB. Figure 6 presents the results, in which the R-QSS-SDF-DPD algorithm performs almost the same as the QSS-SDF-DPD algorithm, but its complexity is much lower (as shown in Figure 4). It means the proposed algorithm is more efficient, and more suitable for practical scenarios.  Figure 7, corroborate that the R-QSS-SDF-DPD algorithm outperforms the SDF-DPD algorithm, whether ULA or NA is deployed. In particular, the R-QSS-SDF-DPD algorithm for NAs is the closest to the CRB, which benefits from the utilization of the QSS features, though a slightly higher complexity has been brought (according to Figure 4). Simulation 4: Different from the last simulation, the function of RMSE versus the number of frames F is adopted in this simulation experiment. Except for the SNR, which is set to be 5 dB and with F varying from 50 to 400, the other parameters are the same as in simulation 3. As presented in Figure 8, it can be found that all algorithms are positively affected by the number of frames, and the R-QSS-SDF-DPD algorithm always outperforms the SDF-DPD algorithm, which is consistent with the results in simulation 3. Simulation 5: Figure 9 evaluates the resolving probability versus the distance between the two sources (∆x). The SNR is assumed to be 0 dB, and the number of frames is F = 200, each frame consists of T F = 600 snapshots; source 1 is located at q 1 = [0, 520]m, and the other one is located at q 2 = [∆x, 520]m, where ∆x varies from 30 m to 150 m. The other parameters are the same as in simulation 3. The result in Figure 9 demonstrates that the resolving probability of the R-QSS-SDF-DPD algorithm for NAs reaches 100% when ∆x is only 50 m, while the SDF-DPD algorithm for NAs reaches 100% when ∆x is 90 m. Obviously, the advantage of the R-QSS-SDF-DPD algorithm, in terms of sources resolution, has been verified.

Conclusions
In this paper, the R-QSS-SDF-DPD algorithm for NAs is proposed. According to the features of QSS, the new noise-eliminated received signal matrix can be obtained by stacking the vectorized form of the signal covariance matrix for each frame and eliminating the unknown noise. Then, the Khatri-Rao subspace method is adopted to reduce the dimension of the cost function. Thereafter, the dimension-reduced matrix is modified, and the unitary transformation method is employed to release the computational burden. It has been verified by some theorical analysis and simulations that the R-QSS-SDF-DPD algorithm owns a lower complexity than the QSS-SDF-DPD algorithm and performs better than the general SDF-DPD algorithm, in terms of the maximum number of identifiable sources, localization accuracy, and source resolutions.

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