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

Abstract

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.

1. 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 likelihood (ML) estimator was established, and multidimensional search was required for the determination of the source position. Though superior performance can be obtained, this algorithm is impractical, in the case of multiple sources, because of its high complexity. To cope up with this, a decoupled algorithm was proposed in [12], and the alternating projection algorithm was adopted in [13]. Besides, the subspace data fusion (SDF) DPD algorithm was proposed in [14], which can handle the problems of multiple sources better than the ML algorithm. It was actually an extension of the multiple signal classification (MUSIC) [15] algorithm. Except for the information regarding TOA and AOA, the information of FDOA can be also considered in the DPD algorithm. In [16], a ML estimator that contained the time difference of arrival (TDOA) and FDOA was designed, and its high computation load was avoided by the particle filter method. Different from the situation of moving receivers in [16], a moving source was considered in [17]. The delay and Doppler information were exploited in [17], and a new multiple particle filter algorithm was proposed to cope with the difficulty of estimating multiple parameters.

The above DPD algorithms were all designed for general Gaussian signals, while some research has demonstrated that an improved accuracy can be achieved if the specific properties of the source signals are considered. In [18,19], the DPD algorithms for the orthogonal frequency division multiplexing (OFDM) signals were proposed, where the ML estimator was exploited, and their superiority over the general algorithms was verified. Besides, some attempts at the DPD algorithms for non-circular (NC) signals have also been made. An improved SDF DPD algorithm for NC signals was proposed in [20], and its complexity was reduced by devising a Newton-type iterative method. Hereafter, some sparse arrays, such as nested array (NA) [21] and coprime array [22], were employed to obtain a larger array aperture, higher degrees of freedom (DOF), and higher accuracy [23,24,25]. In addition to the OFDM and NC signals, the properties of cyclo-stationary signals [26] can also be considered in the DPD algorithm.

Signals in the real world are always nonstationary, but locally stationary, such as speech and audio signals. These types of signals are called quasi-stationary signals (QSS), and they have stable second-order statistical properties within a short period of time and differ from any other frame [27]. To the best of our knowledge, none of the existing literature on DPD algorithms has considered and exploited the features of QSS. In this paper, which was inspired by related research regarding the direction of arrival (DOA) for QSS [27,28,29], the SDF DPD algorithm of QSS (QSS-SDF-DPD) for NAs is derived first. Moreover, the dimension-reduced matrix in the Khatri–Rao subspace method [29] is modified, and the unitary transformation [30,31] is adopted, so that the real-valued QSS-SDF-DPD (R-QSS-SDF-DPD) algorithm with a lower computational burden is proposed. We summarize the following contributions of this paper:

• 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.
• 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.
• 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. $\left|·\right|$ and $‖·‖$ represent the operation of taking the magnitude and Euclidean norm of a vector, respectively. $E\left\{·\right\}$, $vec\left(·\right)$, $diag\left(·\right)$, $blkdiag\left\{·\right\}$, and $Re\left(·\right)$ represent the operation of expectation, vectorization, diagonalization, block diagonalization, and taking the real part, respectively. $\otimes$ and $\odot$ represent the operation of the Kronecker and Khatri–Rao products, respectively. ${ℂ}^{M\times N}$ and ${ℝ}^{M\times N}$ represent the complex number set and real number set with dimension $M\times N$, respectively. The operation of inverse, conjugate, transpose, and conjugate transpose are represented by ${\left(·\right)}^{-1}$, ${\left(·\right)}^{\ast }$, ${\left(·\right)}^T$, and ${\left(·\right)}^H$, respectively. $\arcsin \left(·\right)$ represents the operation of arcsine. $I_n$, $J_n$, $0_N$, and $1_n$ represent the $n\times n$ identity matrix, row-flipped form of the $n\times n$ identity matrix, $n\times 1$ vector with all zeros, and $n\times 1$ vector with all ones.

2. 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=\left[x_n^v,y_n^v\right] \left(n=1,2,\dots ,N\right)$, and the sources are located at $q_k=\left[x_k,y_k\right]$ $\left(k=1,2,\dots ,K\right)$. 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 $\Theta$, given by [21]:

$$\Theta =\left\{d_i=id|i=0,1,\dots ,M_1-1\right\}\cup \left\{d_j=j(M_1+1)d+M_{1}d|j=0,1,\dots ,M_2-1\right\}$$

(1)

where $d$ denotes the unit adjacent spacing.

Assume the $k\mathrm{th}$ source is impinging on the $n\mathrm{th}$ base station from ${\theta }_{n,k}$; then, the received signal vector intercepted by the $n\mathrm{th}$ base station at the $t\mathrm{th} \left(t=1,2,\dots ,T\right)$ sampling time can be presented by [14]:

$$x_n\left(t\right)={\displaystyle \sum _{k=1}^K{a}_n\left({\theta }_{n,k}\right)s_{n,k}\left(t\right)}+n_n\left(t\right)$$

(2)

where $a_n\left({\theta }_{n,k}\right)={\left[e^{-j2\pi d_1\sin {\theta }_{n,k}/\lambda },\dots ,e^{-j2\pi d_m\sin {\theta }_{n,k}/\lambda },\dots ,e^{-j2\pi d_M\sin {\theta }_{n,k}/\lambda }\right]}^T$ denotes the steering vector of ${\theta }_{n,k}=\arcsin \left(\left(x_k-x_n^v\right)/‖q_k-v_n‖\right)$, which is the DOA from the $k\mathrm{th}$ source to the $n\mathrm{th}$ base station, $d_m\in \Theta \left(m=1,2,\dots ,M\right)$, $\lambda$ is the signal wavelength, $s_{n,k}\left(t\right)$ represents the envelope of $k\mathrm{th}$ source incident on the $n\mathrm{th}$ base station, and $n_l\left(t\right)$ is the Gaussian white noise vector.

For the sake of derivation, rewrite Equation (2), in the form of matrix, as [13]:

$$X_n=A_n{S}_n+N_n$$

(3)

where:

$$X_n=\left[x_n\left(1\right),x_n\left(2\right),\dots ,x_n\left(T\right)\right]\in {ℂ}^{M\times T}$$

(4)

$$A_n=\left[a_n\left({\theta }_{n,1}\right),a_n\left({\theta }_{n,2}\right),\dots ,a_n\left({\theta }_{n,K}\right)\right]\in {ℂ}^{M\times K}$$

(5)

$$\begin{array}{cc}{S}_n\hfill & ={\left[s_{n,1},s_{n,2},\cdots ,s_{n,K}\right]}^T\hfill \\ & =\left[\begin{array}{cccc}{s}_{n,1}\left(1\right)& s_{n,1}\left(2\right)& \cdots & s_{n,1}\left(T\right)\\ s_{n,2}\left(1\right)& s_{n,2}\left(2\right)& \cdots & s_{n,2}\left(T\right)\\ ⋮& ⋮& \ddots & ⋮\\ s_{n,K}\left(1\right)& s_{n,K}\left(2\right)& \cdots & s_{n,K}\left(T\right)\end{array}\right]\in {ℂ}^{K\times T}\hfill \end{array}$$

(6)

$$N_n=\left[n_n\left(1\right),n_n\left(2\right),\dots ,n_n\left(T\right)\right]\in {ℂ}^{M\times T}$$

(7)

Considering that the source signals conform to the properties of QSS [36,37,38], we assume $s_{n,k} \left(k=1,2,\dots ,K\right)$ contains $F$ frames of signals, and each frame source signal $s_{n,k,f} \left(f=1,2,\dots ,F\right)$ contains $T_F$ snapshots ($T=F{T}_F$); then, $s_{n,1,f}$ satisfies [27]:

$$E\left\{{\left|s_{n,k,f}\right|}^2\right\}=g_{n,k,f}$$

(8)

Then, the $f\mathrm{th}$ frame received signals of the $n\mathrm{th}$ base station, which can be expressed as:

$$\begin{array}{cc}{X}_{n,f}\hfill & =\left[x_n\left((f-1)T_F+1\right),x_n\left((f-1)T_F+2\right),\dots ,x_n\left(f{T}_f\right)\right]\hfill \\ & =A_n{S}_{n,f}+N_{n,f}\hfill \end{array}$$

(9)

where $S_{n,f}={\left[s_{n,1,f},s_{n,2,f},\dots ,s_{n,K,f}\right]}^T$ denotes the $f\mathrm{th}$ frame source signal matrix, and $N_{n,f}$ is the corresponding additive noise.

According to Equations (8) and (9), a local covariance matrix can be defined by [27]:

$$R_{n,f}=E\left\{X_{n,f}{X}_{n,f}^H\right\}=A_n{\Lambda }_{n,f}{A}_n^H+C_n$$

(10)

where ${\Lambda }_{n,f}=diag\left\{g_{n,1,f},g_{n,2,f},\dots ,g_{n,K,f}\right\}$ denotes the local source covariance matrix of the $f\mathrm{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:

$$y_{n,f}=vec\left(R_{n,f}\right)=vec\left(A_n{\Lambda }_{n,f}{A}_n^H+C_n\right)=\left(A_n^{\ast }\odot A_n\right)g_{n,f}+vec\left(C_n\right)$$

(11)

where $g_{n,f}={\left[g_{n,1,f},g_{n,2,f},\dots ,g_{n,K,f}\right]}^T\in {ℝ}^{K\times 1}$ is a real-value vector consisting of the diagonal elements in ${\Lambda }_{n,f}$.

Interestingly, a new received signal vector $y_{n,f}$ is obtained, whose steering matrix is $A_n^{\ast }\odot A_n$ and noise vector is $vec\left(C_n\right)$. By stacking the vectorization forms of all local covariance matrices, the new received signal matrix $Y_n\in {ℂ}^{M^2\times F}$ is obtained:

$$Y_n=\left[y_{n,1},y_{n,2},\dots ,y_{n,F}\right]=\left(A_n^{\ast }\odot A_n\right)G_n+vec\left(C_n\right)1_F^T$$

(12)

where $G_n=\left[g_{n,1},g_{n,2},\dots ,g_{n,F}\right]\in {ℝ}^{K\times 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.

Define an orthogonal projector matrix $H^{\perp }=I_F-\left(1/F\right)1_F{1}_F^T\in {ℝ}^{F\times F}$, when post-multiplying $Y_n$ by it, the unknown noise in Equation (12) can be eliminated, and the new noise-eliminated received signal matrix ${\tilde{Y}}_n\in {ℂ}^{M^2\times F}$ can be given by [27]:

$$\begin{array}{cc}{\tilde{Y}}_n\hfill & =Y_n{H}^{\perp }=\left(A_n^{\ast }\odot A_n\right)G_n{H}^{\perp }+vec\left(C_n\right)1_F^T{H}^{\perp }\hfill \\ & =\left(A_n^{\ast }\odot A_n\right)G_n{H}^{\perp }+vec\left(C_n\right)1_F^T\left(I_M-\left(1/F\right)1_F{1}_F^T\right)=\left(A_n^{\ast }\odot A_n\right)G_n{H}^{\perp }\hfill \end{array}$$

(13)

According to Assumption A5 in [27], $G_n$ is a full row-rank matrix if $F\ge K+1$, which means $rank\left(G_n^T\right)=K$. Note that $H^{\perp }$ is non-singular, and $G_n{H}^{\perp }$ is a full row-rank matrix. The noise subspace can be obtained after performing singular value decomposition (SVD) on ${\tilde{Y}}_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 ${\tilde{Y}}_n$ is time-consuming. To avoid this problem, the R-QSS-SDF-DPD will be proposed in the next section.

3. The Proposed Algorithm

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

3.1. Khatri–Rao Subspace Method for NA

For a uniform linear array (ULA) with $L$ elements, the new steering matrix $A_n^{\ast }\odot A_n$ can be characterized as [27]:

$$A^{\ast }\odot A=\Gamma B\in {ℂ}^{L^2\times K}$$

(14)

where:

$$\Gamma =\left[\begin{array}{ccccccc}0& \cdots & 0& 1& 0& \cdots & 0\\ 0& \cdots & 0& 0& 1& \cdots & 0\\ ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& \cdots & 0& 0& 0& \cdots & 1\\ 0& \cdots & 1& 0& 0& \cdots & 0\\ 0& \cdots & 0& 1& 0& \cdots & 0\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& \cdots & 0& 0& \cdots & 1& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 1& 0& \cdots & 0& 0& \cdots & 0\\ 0& 1& \cdots & 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 1& 0& \cdots & 0\end{array}\right]\in {ℝ}^{L^2\times \left(2L-1\right)}$$

(15)

$B=\left[b\left({\theta }_1\right),b\left({\theta }_2\right),\dots ,b\left({\theta }_K\right)\right]\in {ℂ}^{\left(2L-1\right)\times K}$ denotes the dimension-reduced virtual array steering matrix, and $b\left({\theta }_k\right)={\left[e^{j2\pi \left(L-1\right)d\sin {\theta }_k/\lambda },\dots ,e^{j2\pi d\sin {\theta }_k/\lambda },1,e^{-j2\pi d\sin {\theta }_k/\lambda },\dots ,e^{-j2\pi \left(L-1\right)d\sin {\theta }_k/\lambda }\right]}^T$ denotes the virtual steering vector of ${\theta }_k$, which is the DOA of the $k\mathrm{th}$ 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_{NA}-1 \left(R_{NA}=M_2\left(M_1+1\right)\right)$ consecutive distinct elements can be obtained after the operation of Khatri–Rao product $A_n^{\ast }\odot 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]:

$$A_n=P{\overline{B}}_n$$

(16)

where $P={\left[p_{d_1},p_{d_2},\dots ,p_{d_M}\right]}^T\in {ℝ}^{M\times R_{NA}}$ is the selection matrix, and $p_{d_m}\in {ℝ}^{R_{NA}\times 1}$ is a vector with 1 at the $j\mathrm{th} \left(j=1+d_m/d\right)$ entry and 0 elsewhere [29]. ${\overline{B}}_n=\left[{\overline{b}}_{n,1},{\overline{b}}_{n,2},\dots ,{\overline{b}}_{n,K}\right]$ $\in {ℂ}^{R_{NA}\times K}$ represents the array steering matrix of the virtual ULA, and ${\overline{b}}_{n,k}={\left[1,e^{-j2\pi d\sin {\theta }_{n,k}/\lambda },\dots ,e^{-j2\pi \left(R_{NA}-1\right)d\sin {\theta }_{n,k}/\lambda }\right]}^T\in {ℂ}^{R_{NA}\times 1}$ represents the corresponding steering vector of the $k\mathrm{th}$ source.

Substitute Equations (16) and (14) into Equation (13); then, ${\tilde{Y}}_n$ can be rewritten as [29]:

$$\begin{array}{cc}{\tilde{Y}}_n\hfill & =\left(A_n^{\ast }\odot A_n\right)G_n{H}^{\perp }=\left({\left(P{\overline{B}}_n\right)}^{\ast }\odot \left(P{\overline{B}}_n\right)\right)G_n{H}^{\perp }\hfill \\ & =\left(P\otimes P\right)\left({\overline{B}}_n^{\ast }\odot {\overline{B}}_n\right)G_n{H}^{\perp }=\left(P\otimes P\right)\left(\tilde{\Gamma }{\tilde{B}}_n\right)G_n{H}^{\perp }\hfill \end{array}$$

(17)

where $\tilde{\Gamma }\in {ℝ}^{R_{NA}{}^2\times \left(2R_{NA}-1\right)}$ is the dimension-reduced matrix with the form of Equation (15), ${\tilde{B}}_n=\left[{\tilde{b}}_n\left({\theta }_{n,1}\right),{\tilde{b}}_n\left({\theta }_{n,2}\right),\dots ,{\tilde{b}}_n\left({\theta }_{n,K}\right)\right]\in {ℂ}^{\left(2R_{NA}-1\right)\times K}$ denotes the dimension-reduced virtual array steering matrix with the form of $B$ in Equation (14), and its $k\mathrm{th}$ column vector ${\tilde{b}}_n\left({\theta }_{n,k}\right)\in$ ${ℂ}^{\left(2R_{NA}-1\right)\times 1}$ represents the corresponding steering vector with the form of $b\left({\theta }_k\right)$ in $B$.

According to the definition of $\tilde{\Gamma }$ and $P$, they are both column orthogonal, so $\left(P\otimes P\right)\tilde{\Gamma }$ is also column orthogonal, which is easy to verify [29]. Define $\tilde{V}={\left(\left(P\otimes P\right)\tilde{\Gamma }\right)}^T\left(P\otimes P\right)\tilde{\Gamma }\in {ℝ}^{\left(2R_{NA}-1\right)\times \left(2R_{NA}-1\right)}$, and the dimension of ${\tilde{Y}}_n$ can be reduced after a linear transformation, which is given by [29]:

$$\begin{array}{cc}{\stackrel{⌣}{Y}}_n& =\tilde{W}{\tilde{Y}}_n={\tilde{V}}^{-1/2}{\left(\left(P\otimes P\right)\tilde{\Gamma }\right)}^T{\tilde{Y}}_n={\tilde{V}}^{-1/2}{\left(\left(P\otimes P\right)\tilde{\Gamma }\right)}^T\left(P\otimes P\right)\left(\tilde{\Gamma }{\tilde{B}}_n\right)G_n{H}^{\perp }\hfill \\ & ={\tilde{V}}^{1/2}{\tilde{B}}_n{G}_n{H}^{\perp }\in {ℂ}^{\left(2R_{NA}-1\right)\times F}\hfill \end{array}$$

(18)

where ${\stackrel{⌣}{Y}}_n$ is the dimension-reduced, noise-eliminated received signal matrix, and $\tilde{W}={\tilde{V}}^{-1/2}{\left(\left(P\otimes P\right)\tilde{\Gamma }\right)}^T$ is the dimension-reduced matrix. As presented in Equation (18), the dimension of the original received signal matrix is reduced from $R_{NA}^2$ to $2R_{NA}-1$.

It can be seen from Equation (18) that a virtual noise-eliminated received signal matrix ${\stackrel{⌣}{Y}}_n$ is obtained, whose array steering matrix is ${\tilde{V}}^{1/2}{\tilde{B}}_n$ and source signal matrix is $G_n{H}^{\perp }\in {ℂ}^{K\times F}$. Therefore, the covariance matrix of the dimension-reduced, noise-eliminated received signal matrix ${\stackrel{⌣}{R}}_n\in {ℂ}^{\left(2R_{NA}-1\right)\times \left(2R_{NA}-1\right)}$ can be calculated by:

$${\widehat{R}}_n={\stackrel{⌣}{Y}}_n{\stackrel{⌣}{Y}}_n^H/F$$

(19)

After the eigenvalue decomposition of ${\widehat{R}}_n$, we obtain the corresponding signal subspace ${\widehat{E}}_n^S$ and noise subspace ${\widehat{E}}_n^N$, which are given by:

$${\widehat{R}}_n={\widehat{E}}_n^S{\widehat{\Lambda }}_n^S{\left({\widehat{E}}_n^S\right)}^H+{\widehat{E}}_n^N{\widehat{\Lambda }}_n^N{\left({\widehat{E}}_n^N\right)}^H$$

(20)

where ${\widehat{\Lambda }}_n^S=diag\left\{{\sigma }_{n,1}^2,{\sigma }_{n,2}^2,\dots ,{\sigma }_{n,K}^2\right\}$ and ${\widehat{\Lambda }}_n^N=diag\left\{{\sigma }_{n,K+1}^2,{\sigma }_{n,K+2}^2,\dots ,{\sigma }_{n,2R_{NA}-1}^2\right\}$ $\left({\sigma }_{n,1}^2>{\sigma }_{n,2}^2,\right.$ $\left.\dots ,{\sigma }_{n,K}^2>{\sigma }_{n,K+1}^2\ge \dots \ge {\sigma }_{n,2R_{NA}-1}^2\right)$ denote the diagonal matrices made up of the maximum $K$ eigenvalues and remaining ones, respectively. The corresponding eigenvectors form ${\widehat{E}}_n^S$ and ${\widehat{E}}_n^N$, respectively.

Hence, according to the SDF algorithm in [14], the spectrum function of the QSS-SDF-DPD algorithm for NAs $s{p}_{QSS-SDF}\left(p\right)$ can be constructed by:

$$s{p}_{QSS-SDF}\left(p\right)=\frac{1}{{\displaystyle \sum _{n=1}^N{\tilde{b}}_n{\left(\theta \left(p\right)\right)}^H{\tilde{V}}^{1/2}{\widehat{E}}_n^N{\left({\widehat{E}}_n^N\right)}^H{\tilde{V}}^{1/2}{\tilde{b}}_n\left(\theta \left(p\right)\right)}}$$

(21)

where $p$ is the position of search grid point in the search area $\Xi$, and ${\tilde{b}}_n\left(\theta \left(p\right)\right)\in {ℂ}^{\left(2R_{NA}-1\right)\times 1}$ represents the steering vector of $p$ for the $n\mathrm{th}$ base station, which has the same form as $b\left({\theta }_k\right)$ in $B$. After finding the maximum $K$ points of Equation (21), all source positions can be estimated.

3.2. Matrix Reconstruction for Complexity Reduction

Note that high complexity cannot be avoided during the peak search process of $s{p}_{SDF}\left(p\right)$, due to the existence of ${\tilde{V}}^{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 ${\stackrel{⌣}{Y}}_n$ by ${\tilde{V}}^{-1/2}$, we can obtain:

$$\begin{array}{cc}{\stackrel{⌣}{Y}}_n^M\hfill & ={\tilde{V}}^{-1/2}{\stackrel{⌣}{Y}}_n={\tilde{B}}_n{G}_n{H}^{\perp }={\tilde{V}}^{-1}\tilde{V}{\tilde{B}}_n{G}_n{H}^{\perp }={\tilde{V}}^{-1}{\left(\left(P\otimes P\right)\tilde{\Gamma }\right)}^T\left(P\otimes P\right)\tilde{\Gamma }{\tilde{B}}_n{G}_n{H}^{\perp }\hfill \\ & ={\tilde{V}}^{-1}{\left(\left(P\otimes P\right)\tilde{\Gamma }\right)}^T{\tilde{Y}}_n={\tilde{W}}^M{\tilde{Y}}_n\hfill \end{array}$$

(22)

where ${\stackrel{⌣}{Y}}_n^M$ is the modified, dimension-reduced, noise-eliminated received signal matrix, and ${\tilde{W}}^M={\tilde{V}}^{-1}{\left(\left(P\otimes P\right)\tilde{\Gamma }\right)}^T\in {ℝ}^{\left(2R_{NA}-1\right)\times M^2}$ is the modified dimension-reduced matrix.

Interestingly and coincidentally, for the NA, ${\tilde{W}}^M$ is equal to the Equation (9) in [39], which is expressed by:

$$\begin{array}{cc}{\left[{\tilde{W}}^M\right]}_{r,j+(k-1)M}\hfill & =\left\{\begin{array}{l}\frac{1}{\omega \left({\Delta }_{j,k}\right)}, {\Delta }_{j,k}=r-R_{NA}\\ 0, \mathrm{otherwise}\end{array}\right.\hfill \\ & r=1,2,\dots ,2R_{NA}-1,j,k=1,2,\dots ,M,\end{array}$$

(23)

where ${\Delta }_{j,k}=\left(d_j-d_k\right)/d \left(d_j,d_k\in \Theta \right)$ represents the position difference between $d_j$ and $d_k$, and $\omega \left({\Delta }_{j,k}\right)$ represents the number of pairs $\left(d_j,d_k\right)$, whose difference $\left(d_j-d_k\right)/d$ is equal to ${\Delta }_{j,k}$.

According to Equation (22), ${\stackrel{⌣}{Y}}_n^M$ can be treated as a modified noise-eliminated received signal matrix, whose array steering matrix is ${\tilde{B}}_n$ and source signal matrix is $G_n{H}^{\perp }\in {ℝ}^{K\times F}$. Note that, compared with the virtual array steering matrix ${\tilde{V}}^{1/2}{\tilde{B}}_n$ before modification, ${\tilde{V}}^{1/2}$ is eliminated, which means the matrix ${\tilde{V}}^{1/2}$ would be removed from $s{p}_{SDF}\left(p\right)$. Thus, it partially reduces the computational burden.

As the virtual array configuration forming ${\tilde{B}}_n$ is centrosymmetric, the unitary transformation method [30,31] can be directly adopted.

Define the unitary matrix by [30]:

$$U_{2R_{NA}-1}=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_{R_{NA}-1}& 0_{R_{NA}-1}& j{J}_{R_{NA}-1}\\ 0_{R_{NA}-1}^T& \sqrt{2}& 0_{R_{NA}-1}^T\\ J_{R_{NA}-1}& 0_{R_{NA}-1}& -j{I}_{R_{NA}-1}\end{array}\right]$$

(24)

Then, a real-valued matrix ${\stackrel{⌣}{Y}}_n^{MR}\in {ℝ}^{\left(2R_{NA}-1\right)\times F}$ can be obtained by the unitary transformation given by:

$${\stackrel{⌣}{Y}}_n^{MR}=U_{2R_{NA}-1}^H{\stackrel{⌣}{Y}}_n^M$$

(25)

Correspondingly, a real-valued covariance matrix ${\widehat{R}}_n^{MR}\in {ℝ}^{\left(2R_{NA}-1\right)\times \left(2R_{NA}-1\right)}$ can be expressed by:

$${\widehat{R}}_n^{MR}={\stackrel{⌣}{Y}}_n^{MR}{\left({\stackrel{⌣}{Y}}_n^{MR}\right)}^H/F$$

(26)

Hence, after performing eigenvalue decomposition on ${\widehat{R}}_n^{MR}$, we obtain:

$${\widehat{R}}_n^{MR}={\widehat{E}}_n^{MRS}{\widehat{\Lambda }}_n^{RS}{\left({\widehat{E}}_n^{MRS}\right)}^H+{\widehat{E}}_n^{MRN}{\widehat{\Lambda }}_n^{RN}{\left({\widehat{E}}_n^{MRN}\right)}^H$$

(27)

where ${\widehat{E}}_n^{MRS}$ is the real-valued signal subspace made up of eigenvectors corresponding to the largest $K$ eigenvalues, ${\widehat{E}}_n^{MRN}$ is the real-valued noise subspace made up of the remaining eigenvectors, and ${\left({\sigma }_{n,1}^R\right)}^2>{\left({\sigma }_{n,2}^R\right)}^2,\dots ,{\left({\sigma }_{n,K}^R\right)}^2>{\left({\sigma }_{n,K+1}^R\right)}^2\ge \dots \ge {\left({\sigma }_{n,2R_{NA}-1}^R\right)}^2$ are eigenvalues of ${\widehat{R}}_n^{MR}$; the first $K$ ones form the diagonal matrix ${\widehat{\Lambda }}_n^{RS}$, and the last $2R_{NA}-K-1$ ones form the diagonal matrix ${\widehat{\Lambda }}_n^{RN}$.

Finally, fuse all the real-value noise subspace of all base stations, and the spectrum function of R-QSS-SDF-DPD algorithm for NAs $s{p}_{R-QSS-SDF}\left(p\right)$ can be expressed by:

$$s{p}_{R-QSS-SDF}\left(p\right)=\frac{1}{{\displaystyle \sum _{n=1}^N{\tilde{b}}_n^R{\left(\theta \left(p\right)\right)}^H{\widehat{E}}_n^{MRN}{\left({\widehat{E}}_n^{MRN}\right)}^H{\tilde{b}}_n^R\left(\theta \left(p\right)\right)}}$$

(28)

where ${\tilde{b}}_n^R\left(\theta \left(p\right)\right)=U_{2R_{NA}-1}^H{\tilde{b}}_n\left(\theta \left(p\right)\right)=\sqrt{2}\left[\cos \left(2\pi \left(R_{NA}-1\right)d\sin \theta \left(p\right)/\lambda \right),\dots ,\right.$ $\cos \left(2\pi d\sin \theta \left(p\right)/\lambda \right),1/\sqrt{2}\left.,\sin \left(2\pi d\sin \theta \left(p\right)/\lambda \right),\dots ,\sin \left(2\pi \left(R_{NA}-1\right)d\sin \theta \left(p\right)/\lambda \right)\right]\in {ℝ}^{\left(2R_{NA}-1\right)\times 1}$ is the real-valued peak search steering vector. Then, find the maximum $K$ points of $s{p}_{MR-QSS-SDF}\left(p\right)$, and take the corresponding coordinate positions as the estimates of all source positions.

3.3. 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 3: Eliminate the unknown noise, according to Equation (13).

Step 4: Obtain the original dimension-reduced, noise-eliminated received signal matrix ${\stackrel{⌣}{Y}}_n$ by Equation (18).

Step 5: Calculate the covariance matrix of ${\stackrel{⌣}{Y}}_n$ by Equation (19), and obtain the noise subspace matrix ${\widehat{E}}_n^N$ by Equation (20).

Step 6: Construct the original spectrum function of QSS-SDF-DPD $s{p}_{QSS-SDF}\left(p\right)$ 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 ${\stackrel{⌣}{Y}}_n^M$ according to Equation (22).

Step 8: Obtain the real-valued matrix ${\stackrel{⌣}{Y}}_n^{MR}$, according to Equation (25).

Step 9: Calculate the real-valued covariance matrix ${\widehat{R}}_n^{MR}$ by Equation (26), and obtain the real-valued noise subspace matrix ${\widehat{E}}_n^{MRN}$ by Equation (27).

Step 10: Construct the spectrum function of R-QSS-SDF-DPD $s{p}_{R-QSS-SDF}\left(p\right)$ by Equation (28), and estimate all source positions by finding the maximum $K$ points.

4. Performance Analysis

4.1. CRB

Based on the derivation results in [24,40], the CRB of DPD for NAs can be expressed by:

$$CRB\left(q\right)=\frac{{\sigma }^2}{2}{\left\{{\displaystyle \sum _{t=1}^{T}Re\left[S_t^H{D}^H{\Pi }_A^{\perp }D{S}_t\right]}\right\}}^{-1}$$

(29)

where ${\sigma }^2$ is the variance of the noise, and:

$$q={\left[q_1,q_2,\dots ,q_K\right]}^T$$

(30)

$$S_t=I_{2K}\otimes s_t$$

(31)

$$s_t={\left[s_{1,t}^T,s_{2,t}^T,\dots ,s_{N,t}^T\right]}^T$$

(32)

$${\Pi }_A^{\perp }=I-A{\left(A^{H}A\right)}^{-1}{A}^H$$

(33)

$$A=blkdiag\left\{A_1,A_2,\dots ,A_N\right\}$$

(34)

$$D=\left[\frac{\partial A}{\partial x_1},\frac{\partial A}{\partial y_1},\frac{\partial A}{\partial x_2},\frac{\partial A}{\partial y_2},\dots ,\frac{\partial A}{\partial x_K},\frac{\partial A}{\partial y_K}\right]$$

(35)

Where $\frac{\partial A}{\partial x_k}$ and $\frac{\partial A}{\partial y_k}$ denote the partial derivatives of $A$, with respect to $x_k$ and $y_k$, respectively.

4.2. 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. Complexity of algorithms.

Algorithm	Major Steps	Complexity (Real-Valued Multiplication Operation)
The SDF-DPD	Obtain covariance matrix	$O\left\{4N{M}^{2}T\right\}$
	Spatial smoothing	$O\left\{4N{R}_{NA}^3\right\}$
	EVD	$O\left\{4N{R}_{NA}^3\right\}$
	Spectral peak search	$O\{4N(R_{NA}^2(R_{NA}-K)+N_x{N}_y{R}_{NA}(R_{NA}+1))\}$
	Total	$O\{4N(M^{2}T+3R_{NA}^3-K{R}_{NA}^2+N_x{N}_y{R}_{NA}(R_{NA}+1))\}$
The proposed QSS-SDF-DPD	Obtain $R_{n,f}$	$O\left\{4N{M}^{2}T\right\}$
	Obtain ${\tilde{Y}}_n$	$O\left\{2N{M}^2{F}^2\right\}$
	Obtain ${\stackrel{⌣}{Y}}_n$	$O\{2N{M}^{2}F(2R_{NA}-1)\}$
	Obtain ${\widehat{R}}_n$	$O\left\{4NF{(2R_{NA}-1)}^2\right\}$
	EVD	$O\left\{4N{(2R_{NA}-1)}^3\right\}$
	Spectral peak search	$O\{2N(3{(2R_{NA}-1)}^2(2R_{NA}-K-1)+4N_x{N}_y{R}_{NA}(2R_{NA}-1))\}$
	Total	$O\{2N(2M^{2}T+M^2{F}^2+M^{2}F(2R_{NA}-1)+(2F-3K)$${(2R_{NA}-1)}^2+5{(2R_{NA}-1)}^3+4N_x{N}_y{R}_{NA}(2R_{NA}-1))\}$
The proposed R-QSS-SDF-DPD	Obtain $R_{n,f}$	$O\left\{4N{M}^{2}T\right\}$
	Obtain ${\tilde{Y}}_n$	$O\left\{2N{M}^2{F}^2\right\}$
	Obtain ${\stackrel{⌣}{Y}}_n^M$	$O\{2N{M}^{2}F(2R_{NA}-1)\}$
	Obtain ${\stackrel{⌣}{Y}}_n^{MR}$	$O\left\{4NF{\left(2R_{NA}-1\right)}^2\right\}$
	Obtain ${\widehat{R}}_n^{MR}$	$O\left\{NF{(2R_{NA}-1)}^2\right\}$
	EVD	$O\left\{N{(2R_{NA}-1)}^3\right\}$
	Spectral peak search	$O\{N({(2R_{NA}-1)}^2(2R_{NA}-K-1)+2N_x{N}_y{R}_{NA}(2R_{NA}-1))\}$
	Total	$O\{N(4M^{2}T+2M^2{F}^2+2M^{2}F(2R_{NA}-1)+(5F-K)$${(2R_{NA}-1)}^2+2{(2R_{NA}-1)}^3+2N_x{N}_y{R}_{NA}(2R_{NA}-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\times N_y$ grids, $M=M_1+M_2$, and $R_{NA}=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.

4.3. 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.

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

5. 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:

$$\mathrm{RMSE}=\frac{1}{K}{\displaystyle \sum _{k=1}^K\sqrt{\frac{1}{N_E}{\displaystyle \sum _{ne=1}^{N_E}\left[{\left({\widehat{x}}_{k,ne}-x_k\right)}^2+{\left({\widehat{y}}_{k,ne}-y_k\right)}^2\right]}}}$$

(36)

where $N_E$ denotes the number of Monte Carlo simulation experiments, and $\left[{\widehat{x}}_{k,ne},{\widehat{y}}_{k,ne}\right]$ denotes the estimate of $\left[x_k,y_k\right]$ in the $ne\mathrm{th}$ 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:

$$P_r=\frac{N_s}{N_E}\times 100\%$$

(37)

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=\left[x_1,y_1\right]m$, while the abscissa of the other $q_2=\left[x_1+\Delta x,y_2\right]m$ is adjusted, where $\Delta x$ denotes the distance of the two sources. For each simulation experiment, if the locations of two sources are estimated and satisfy $\left|{\widehat{x}}_1-x_1\right|<\Delta x/2, \left|{\widehat{x}}_2-x_2\right|<\Delta x/2$ (${\widehat{x}}_1$ and ${\widehat{x}}_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=[-800,400]m$, $q_2=[-200,-600]m$, $q_3=[-400,1000]m$, $q_4=[700,1300]m$, $q_5=[1200,0]m$, $q_6=[1600,400]m$, and $q_7=[800,500]m$; the number of base stations is $N=4$, and they are located at $v_1=[0,2000]m$, $v_2=[-700,-1600]m$, $v_3=[1100,-1200]m$, and $v_4=[2000,2100]m$. The structure of NA is $M_1=2$, $M_2=2$, SNR is 5 dB, and the number of frames is $F=200$; each frame consists of $T_F=600$ snapshots, and the number of Monte Carlo simulation experiments is $500$. As demonstrated in Figure 5, the proposed algorithm can estimate all source positions, even when $K=7>R_{NA}-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.

Simulation 2: Figure 6 depicts the RMSE performance of algorithms versus SNR, where the R-QSS-SDF-DPD algorithm is compared to its version before matrix reconstruction (QSS-SDF-DPD). Here, we consider two conditions, where ULA and NA are deployed. The number of sources is set to $K=2$, $q_1=\left[0,390\right]m$, and $q_2=[250,560]m$; the number of base stations is $N=3$, $v_1=[-700,-300]m$, $v_2=[200,-500]m$, and $v_3=[600,-200]m$. 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.

Simulation 3: Figure 7 depicts the RMSE performance of algorithms versus SNR, where the R-QSS-SDF-DPD algorithm is compared with the SDF-DPD algorithm which ignores the features of QSS. In this simulation experiment, the number of sources is $K=2$, $q_1=[0,390]m$, and $q_2=[60,580]m$; the number of base stations is $N=3$, $v_1=[-700,-300]m$, $v_2=[200,-500]m$, and $v_3=[600,-200]m$. The structure of the array is $M=4$ for ULA and $M_1=2$, $M_2=2$ for NA; the number of frames is $F=400$, each frame consists of $T_F=300$ snapshots, and the SNR varies from −3 dB to 15 dB. The results, depicted in 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 ($\Delta 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=[\Delta x,520]m$, where $\Delta 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 $\Delta x$ is only 50 m, while the SDF-DPD algorithm for NAs reaches 100% when $\Delta x$ is 90 m. Obviously, the advantage of the R-QSS-SDF-DPD algorithm, in terms of sources resolution, has been verified.

6. 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.