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.

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. and represent the operation of taking the magnitude and Euclidean norm of a vector, respectively. , , , , and 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. and represent the complex number set and real number set with dimension , respectively. The operation of inverse, conjugate, transpose, and conjugate transpose are represented by , , , and , respectively. represents the operation of arcsine. , , , and represent the identity matrix, row-flipped form of the identity matrix, vector with all zeros, and vector with all ones.

2. System Model

Consider the two-dimension scenario presented in Figure 1, where the (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 base stations, which are equipped with a NA. As the location of the base stations are known, assume they are located at , and the sources are located at . The specific structure of the element NA that is exploited in this scenario is shown in Figure 2, where the first and second levels consist of and elements (), respectively. The place of all physical array elements can be included in a set , given by [21]: where denotes the unit adjacent spacing.

Figure 1

The localization scenario for multiple sources with multiple NAs.

Figure 2

The structure of NA.

Assume the source is impinging on the base station from ; then, the received signal vector intercepted by the base station at the sampling time can be presented by [14]: where denotes the steering vector of , which is the DOA from the source to the base station, , is the signal wavelength, represents the envelope of source incident on the base station, and is the Gaussian white noise vector.

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

Considering that the source signals conform to the properties of QSS [36,37,38], we assume contains frames of signals, and each frame source signal contains snapshots (); then, satisfies [27]:

Then, the frame received signals of the base station, which can be expressed as: where denotes the frame source signal matrix, and is the corresponding additive noise.

According to Equations (8) and (9), a local covariance matrix can be defined by [27]: where denotes the local source covariance matrix of the frame, and is the spatial noise covariance.

According to Property 1 in [27], the vectorization of in Equation (10) can be expressed as: where is a real-value vector consisting of the diagonal elements in .

Interestingly, a new received signal vector is obtained, whose steering matrix is and noise vector is . By stacking the vectorization forms of all local covariance matrices, the new received signal matrix is obtained: where is the new signal matrix.

Based on Assumption A4 in [27], the covariance of the source signal is time-varying. Therefore, can be treated as a new incoherent source signal matrix with 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.

Figure 3

The DOF of the general signal model and QSS model.

Define an orthogonal projector matrix , when post-multiplying by it, the unknown noise in Equation (12) can be eliminated, and the new noise-eliminated received signal matrix can be given by [27]:

According to Assumption A5 in [27], is a full row-rank matrix if , which means . Note that is non-singular, and is a full row-rank matrix. The noise subspace can be obtained after performing singular value decomposition (SVD) on 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 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 elements, the new steering matrix can be characterized as [27]: where: denotes the dimension-reduced virtual array steering matrix, and denotes the virtual steering vector of , which is the DOA of the 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 .

For the two-level NA exploited in this paper, a total of consecutive distinct elements can be obtained after the operation of Khatri–Rao product . 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 is the selection matrix, and is a vector with 1 at the entry and 0 elsewhere [29]. represents the array steering matrix of the virtual ULA, and represents the corresponding steering vector of the source.

Substitute Equations (16) and (14) into Equation (13); then, can be rewritten as [29]: where is the dimension-reduced matrix with the form of Equation (15), denotes the dimension-reduced virtual array steering matrix with the form of in Equation (14), and its column vector represents the corresponding steering vector with the form of in .

According to the definition of and , they are both column orthogonal, so is also column orthogonal, which is easy to verify [29]. Define , and the dimension of can be reduced after a linear transformation, which is given by [29]: where is the dimension-reduced, noise-eliminated received signal matrix, and is the dimension-reduced matrix. As presented in Equation (18), the dimension of the original received signal matrix is reduced from to .

It can be seen from Equation (18) that a virtual noise-eliminated received signal matrix is obtained, whose array steering matrix is and source signal matrix is . Therefore, the covariance matrix of the dimension-reduced, noise-eliminated received signal matrix can be calculated by:

After the eigenvalue decomposition of , we obtain the corresponding signal subspace and noise subspace , which are given by: where and denote the diagonal matrices made up of the maximum eigenvalues and remaining ones, respectively. The corresponding eigenvectors form and , respectively.

Hence, according to the SDF algorithm in [14], the spectrum function of the QSS-SDF-DPD algorithm for NAs can be constructed by: where is the position of search grid point in the search area , and represents the steering vector of for the base station, which has the same form as in . After finding the maximum 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 , due to the existence of 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 by , we can obtain: where is the modified, dimension-reduced, noise-eliminated received signal matrix, and is the modified dimension-reduced matrix.

Interestingly and coincidentally, for the NA, is equal to the Equation (9) in [39], which is expressed by: where represents the position difference between and , and represents the number of pairs , whose difference is equal to .

According to Equation (22), can be treated as a modified noise-eliminated received signal matrix, whose array steering matrix is and source signal matrix is . Note that, compared with the virtual array steering matrix before modification, is eliminated, which means the matrix would be removed from . Thus, it partially reduces the computational burden.

As the virtual array configuration forming is centrosymmetric, the unitary transformation method [30,31] can be directly adopted.

Define the unitary matrix by [30]:

Then, a real-valued matrix can be obtained by the unitary transformation given by:

Correspondingly, a real-valued covariance matrix can be expressed by:

Hence, after performing eigenvalue decomposition on , we obtain: where is the real-valued signal subspace made up of eigenvectors corresponding to the largest eigenvalues, is the real-valued noise subspace made up of the remaining eigenvectors, and are eigenvalues of ; the first ones form the diagonal matrix , and the last ones form the diagonal matrix .

Finally, fuse all the real-value noise subspace of all base stations, and the spectrum function of R-QSS-SDF-DPD algorithm for NAs can be expressed by: where   is the real-valued peak search steering vector. Then, find the maximum points of , 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 by Equation (10).

Step 2: Construct the vectorization forms of all , and stack them together to obtain the new received signal matrix 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 by Equation (18).

Step 5: Calculate the covariance matrix of by Equation (19), and obtain the noise subspace matrix by Equation (20).

Step 6: Construct the original spectrum function of QSS-SDF-DPD by Equation (21), and estimate all source positions by finding the maximum points.

Step 7: Modify the dimension-reduced matrix to obtain the modified noise-eliminated received signal matrix according to Equation (22).

Step 8: Obtain the real-valued matrix , according to Equation (25).

Step 9: Calculate the real-valued covariance matrix by Equation (26), and obtain the real-valued noise subspace matrix by Equation (27).

Step 10: Construct the spectrum function of R-QSS-SDF-DPD by Equation (28), and estimate all source positions by finding the maximum 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: where is the variance of the noise, and: Where and denote the partial derivatives of , with respect to and , 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 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 grids, , and .

Table 1

Complexity of algorithms.

Algorithm	Major Steps	Complexity (Real-Valued Multiplication Operation)
The SDF-DPD	Obtain covariance matrix	O{4NM^2T}
	Spatial smoothing	O{4NRNA3}
	EVD	O{4NRNA3}
	Spectral peak search	O{4N(RNA2(R_NA−K)+N_xN_yR_NA(R_NA+1))}
	Total	O{4N(M^2T+3RNA3−KRNA2+N_xN_yR_NA(R_NA+1))}
The proposed QSS-SDF-DPD	Obtain R_n,f	O{4NM^2T}
	Obtain Y˜_n	O{2NM^2F^2}
	Obtain Y⌣_n	O{2NM^2F(2R_NA−1)}
	Obtain R^_n	O{4NF(2R_NA−1)^2}
	EVD	O{4N(2R_NA−1)^3}
	Spectral peak search	O{2N(3(2R_NA−1)^2(2R_NA−K−1)+4N_xN_yR_NA(2R_NA−1))}
	Total	O{2N(2M^2T+M^2F^2+M^2F(2R_NA−1)+(2F−3K) (2R_NA−1)^2+5(2R_NA−1)^3+4N_xN_yR_NA(2R_NA−1))}
The proposed R-QSS-SDF-DPD	Obtain R_n,f	O{4NM^2T}
	Obtain Y˜_n	O{2NM^2F^2}
	Obtain Y⌣nM	O{2NM^2F(2R_NA−1)}
	Obtain Y⌣nMR	O{4NF2R_NA−1^2}
	Obtain R^nMR	O{NF(2R_NA−1)^2}
	EVD	O{N(2R_NA−1)^3}
	Spectral peak search	O{N((2R_NA−1)^2(2R_NA−K−1)+2N_xN_yR_NA(2R_NA−1))}
	Total	O{N(4M^2T+2M^2F^2+2M^2F(2R_NA−1)+(5F−K) (2R_NA−1)^2+2(2R_NA−1)^3+2N_xN_yR_NA(2R_NA−1))}

Figure 4

Comparison of complexity.

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 ( is the actual number of sources, and 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: where denotes the number of Monte Carlo simulation experiments, and denotes the estimate of in the simulation experiment. 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 [41] is defined by: where 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 , while the abscissa of the other is adjusted, where denotes the distance of the two sources. For each simulation experiment, if the locations of two sources are estimated and satisfy ( and 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 , and they are located at , , , , , , and ; the number of base stations is , and they are located at , , , and . The structure of NA is , , SNR is 5 dB, and the number of frames is ; each frame consists of snapshots, and the number of Monte Carlo simulation experiments is . As demonstrated in Figure 5, the proposed algorithm can estimate all source positions, even when , 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.

Figure 5

Scatter plots of the R-QSS-SDF-DPD algorithm for NAs.

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 , , and ; the number of base stations is , , , and . The structure of array is for ULA, and , for NA. The number of frames is , each frame consists of 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 6

Comparison of QSS-SDF-DPD and R-QSS-SDF-DPD algorithms by RMSE versus SNR.

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 , , and ; the number of base stations is , , , and . The structure of the array is for ULA and , for NA; the number of frames is , each frame consists of 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).

Figure 7

Comparison of SDF-DPD and R-QSS-SDF-DPD algorithms by RMSE versus SNR.

Simulation 4: Different from the last simulation, the function of RMSE versus the number of frames is adopted in this simulation experiment. Except for the SNR, which is set to be 5 dB and with 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.

Figure 8

Comparison of SDF-DPD and R-QSS-SDF-DPD algorithms by RMSE versus the number of frames.

Simulation 5: Figure 9 evaluates the resolving probability versus the distance between the two sources (). The SNR is assumed to be 0 dB, and the number of frames is , each frame consists of snapshots; source 1 is located at , and the other one is located at , where 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 is only 50 m, while the SDF-DPD algorithm for NAs reaches 100% when is 90 m. Obviously, the advantage of the R-QSS-SDF-DPD algorithm, in terms of sources resolution, has been verified.

Figure 9

Comparison of SDF-DPD and R-QSS-SDF-DPD algorithms by resolving probability versus .

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.