An Improved Trilinear Model-Based Angle Estimation Method for Co-Prime Planar Arrays

Angle estimation methods in two-dimensional co-prime planar arrays have been discussed mainly based on peak searching and sparse recovery. Peak searching methods suffer from heavy computational complexity and sparse recovery methods face some problems in selecting the regularization parameters. In this paper, we propose an improved trilinear model-based method for angle estimation for co-prime planar arrays in the view of trilinear decomposition, namely parallel factor analysis. Due to the principle of trilinear decomposition, our method does not require peak searching and can conduct auto-pairing easily, which can reduce the computational loads and avoid parameter selection problems. Furthermore, we exploit the virtual array concept of the whole co-prime planar array through the cross-correlation matrix obtained from the received signal data and present a matrix reconstruction method using the Khatri–Rao product to tackle the matrix rank deficiency problem in the virtual array condition. The simulation results show that our proposed method can not only achieve high estimation accuracy with low complexity compared to other similar approaches, but also utilize limited sensor number to implement the angle estimation tasks.


Introduction
Angle estimation for planar arrays, i.e., two-dimensional direction of arrival (2D DOA) estimation, plays an important role in various applications, such as radar, sonar, and wireless communications [1]. For instance, in monitoring the electromagnetic environment application, radar signal parameters are detected in the situations, which are always noisy or changeable in relation to the weather conditions. Therefore, the precise and fast measurement of those parameters with a high accuracy is very important. Some modern technologies of planar arrays based on Artificial Neural Networks [2,3] and different radar localization methods [4] have been proposed. The key problem is to estimate the azimuth and elevation angles correctly and rapidly without angle mismatches. Scholars and researchers have done much work on this problem from different views. The author of Ref. [5] proposed a new computationally efficient cross-correlation based 2D DOA estimation (CODE) method without eigen decomposition for 2-D DOA estimation of non-coherent signals impinging on the L-shaped sensor array. The conjugate symmetry property of the array manifold matrix is utilized in Ref. [6] to increase the effective array aperture, which improves the angle estimation performance. The researchers in Ref. [7] proposed a method based on the matrix reconstruction obtained from the cross-correlation matrix between two sub-arrays, hence reducing the computational complexity. These methods are designed for L-shaped array structures. Studies about 2D DOA estimation methods for uniform rectangular arrays (URAs) also attract significant attention. For example, the author of Ref. [8] introduced a preprocessing transformation matrix and performed real-valued calculations. In Ref. [9], the scholars presented a linear Least Squares method for solving the joint 2D DOA estimation problem under a single source case.
All of the methods mentioned above primarily deal with arrays with uniform spacing equal to half a wavelength of sensor elements, which limits the maximum source number that the array structure can detect. Sparse arrays open up the field of sensor array processing because of the higher degrees of freedom (DOFs) offered by the co-array domain. Recently, a new class of array geometry named co-prime array [10,11] was proposed to increase the DOFs and enhance the performance of the array structure. The DOA methods were presented after obtaining the greater DOFs. There has been a proliferation of studies considering DOA algorithms. The author of Refs. [12,13] introduced spatial smoothing techniques and Toeplitz matrix reconstruction based on the Multiple Signal Classification (MUSIC) algorithm [14] or any other grid-less algorithm [15]. However, these methods only utilize the consecutive parts of the difference co-array and the DOFs are not fully exploited. To overcome the problem of aperture loss, the compressive sensing (CS)-based methods [16][17][18] were developed to improve the estimation resolution. Other works, such as Refs. [19][20][21] processed the signals from a generalized co-prime array, fourth-order difference co-array concept and super resolution theory, respectively. It is regrettable that those methods [12][13][14][15][16][17][18][19][20][21] were developed for co-prime linear arrays (CPLAs) and cannot be exploited in the2D DOA estimation area directly. According to the array structure in the 2D case, there exists a kind of sparse planar array named the co-prime planar array (CPPA). The authors of Refs. [22][23][24] discussed the array structure and DOA estimation methods in 2D conditions, where a lattice generator and 2D spatial smoothing were utilized. However, the method requires spatial smoothing matrix construction based on the peak searching algorithm, such as MUSIC, which leads to a heavy computational load. The author of Ref. [25] employed a partial spectrum peak searching method that is computationally efficient. Nonetheless, they did not exploit the virtual array concept of the CPPA and ignored the aperture extensions. Other methods in Refs. [26,27] also suffers from great computational complexity caused by eigenvalue decomposition or 2D grids searching.
In this paper, we concentrate on the virtual array concept of the CPPA derived from the cross-correlation matrix and construct an improved trilinear model for the angle estimation problem. Rather than estimating the azimuth and elevation angles from decomposing the CPPA into two sub-arrays, we calculate the cross-correlation matrix between the two sub-arrays and acquire the virtual array of the CPPA. Larger DOFs can be achieved by selecting the consecutive parts of the whole virtual array. Moreover, we present a matrix reconstruction method to resolve the problem of matrix deficiency under the virtual array concept. It is remarkable, that based on conventional trilinear decomposition an improved trilinear model is set up by applying a compressed matrix. Detailed analysis and simulation results show that our proposed method yields low computational complexity and better estimation performance. More importantly, our proposed method does not require peak searching or 2D grids sampling, as compared to the aforementioned methods, and can achieve angle auto-pairing.
The rest of the paper is organized as follows Section 2 discusses the array structure and signal model used in this paper. The proposed method is given in Section 3. Simulation and analysis are performed in Section 4. The conclusion is presented in Section 5.

Array Structure and Signal Model
We assume that the geometry of the CPPA is constructed by the following steps: firstly, two uniform linear arrays (ULAs), in which each sub-array contains N sensors having inter-element spacing Md d = λ 2 and 2M sensors having inter-element spacing Nd d = λ 2 , respectively, are concatenated together with the 0th element aligned into one non-uniform linear array (NLA). Suppose that N and M are co-prime integers (M < N); the NLA is actually a kind of structure named CPLA, as depicted in Figure 1a  Secondly, the two CPLAs generated in the first step are displaced along the x direction and the y direction in a Cartesian Rectangular Coordinate System, where the 0 th element of each subarray are put together in the original position. Finally, we extend the arranged CPLAs into the geometry of the CPPA, as shown in Figure 2. From Figure 2 we notice that the CPPA consists of two uniform rectangular arrays (URAs) with sizes of 2 2 M M  and N N  , respectively. The total number of the whole array is  signal source, respectively. In order to exploit the co-prime property and analyze Secondly, the two CPLAs generated in the first step are displaced along the x direction and the y direction in a Cartesian Rectangular Coordinate System, where the 0th element of each subarray are put together in the original position. Finally, we extend the arranged CPLAs into the geometry of the CPPA, as shown in Figure 2. From Figure 2 we notice that the CPPA consists of two uniform rectangular arrays (URAs) with sizes of 2M × 2M and N × N, respectively. The total number of the whole array is   The sensors in the CPPA model are assumed to be identical, omnidirectional, and isotropic. K uncorrelated narrowband far-field signals {s k (t), k = 1, 2, . . . , K} impinge on the CPPA from distinct directions. θ k (k = 1, 2, . . . , K) and φ k (k = 1, 2, . . . , K) represent the elevation angle and the azimuth angle of the kth signal source, respectively. In order to exploit the co-prime property and analyze the DOA problem further, we start our signal model from the array-divided view, utilizing the sub-array from the whole CPPA. We extract the URA with a size of 2M × 2M from the structure and the received lth snapshot data can be expressed as: where A 1 (θ k , φ k ) ∈ C 2M×2M is the manifold matrix for the selected URA and • represents the outer product of two vectors. a 1,x (θ k , φ k ) ∈ C 2M×1 and a 1,y (θ k , φ k ) ∈ C 2M×1 can be denoted as: Here we use the definitions as follows: Then we can rewrite Equations (2) and (3) as: where n 1 (l) ∼ CN 0, σ 2 n I represents the complex-valued additive white Gaussian noise term and σ 2 n indicates the noise power. The vectorization of a matrix will give simplicity in further analysis and discussion, so we can obtain the substituted with Equation (1) as follows: where A 1,y = a 1,y (θ 1 , φ 1 ), a 1,y (θ 2 , φ 2 ), . . . , a 1,y (θ K , φ K ) = a 1,y r 1,y , a 1,y r 2,y , . . . , a 1,y r K,y ∈ includes the signal waveforms under the lth snapshot condition and N 1 (l) = vec(n 1 (l)) ∈ C 4M 2 ×1 . According to the derivations above, we can naturally obtain the received signal expression for another sub-array isolated from the CPPA with a size of N × N as follows: where A 2,y = a 2,y (θ 1 , φ 1 ), a 2,y (θ 2 , φ 2 ), . . . , a 2,y (θ K , φ K ) = a 2,y r 1,y , a 2,y r 2,y , . . . , a 2,y r K,y ∈ is the complex-valued additive white Gaussian noise term with zero-mean and variance σ 2 n .

The Proposed Angle Estimation Method
Unlike those methods exploiting two sub-arrays for peak searching or using predefined sampling 2D grids in DOA estimation, we propose a new method based on an improved trilinear model for the whole CPPA. Firstly, we extend the aperture of the CPPA by utilizing the virtual array concept, which increases the degrees of freedom (DOFs). During the process of aperture extension, we transform the original received signals of two sub-arrays into an equivalent second-order received signal of the virtual URA array; thus the DOA estimation problem for the CPPA can be simplified into the same problem for the virtual URA. Secondly, we construct an improved trilinear model for the equivalent received signal and propose a corresponding angle estimation method through detailed mathematical deduction. The proposed method will enhance the estimation performance by reducing the dimensionality of the conventional trilinear model into the improved one.

Virtual Array Concept For CPPA
As mentioned in Ref. [13], transforming the received data into their second-order (or high-order) statistics is the essence of sparse array signal processing. The theoretical cross-correlation matrix between the received signals of two subarrays decomposed from the CPPA can be calculated as: where ∈ C 4M 2 ×N 2 4M 2 > N 2 (note that due to the co-primness between integers M and N, 4M 2 = N 2 ). Then we can obtain the vectorization of R 12 : where and p = σ 2 1 , σ 2 2 , . . . , σ 2 K T , q = vec(N 12 ). In Equation (12), B behaves like the manifold of a virtual array whose sensor locations in the xoy plane are elaborated by the distinct elements in the set , then it has been proved in Ref. [19] that set S contains all the integers in the range [−(N − 1), MN + M − 1]. Thus, the virtual array sensor will be displaced at a total of (MN + M + N − 1) × (MN + M + N − 1) = (MN + M + N − 1) 2 numbers and the DOFs will be increased dramatically. Figure 3 shows the virtual array concept for the CPPA in view of the equivalent second-order received signal.
, then it has been proved in Ref. [19] that set  contains all the integers in the range Thus, the virtual array sensor will be displaced at a total of      numbers and the DOFs will be increased dramatically. Figure 3 shows the virtual array concept for the CPPA in view of the equivalent second-order received signal.

Matrix Reconstruction Method
According to the analysis of the aforementioned equivalent received signal model, we can whose elements correspond to the continuous rows of B as a virtual array manifold matrix from the location can model the equivalent received signal as follows: where: .. .
and q  represents the deterministic noise vector corresponding to the terms with the same rows as B  in q .

Remarks:
The equivalent received signal exhibits a form similar to that of the true received signal expression. In spite of that, subspace-based DOA methods such as MUSIC or ESPRIT cannot be used on the equivalent

Matrix Reconstruction Method
According to the analysis of the aforementioned equivalent received signal model, we can construct a new matrix B with a size of (MN + M + N − 1) 2 × K whose elements correspond to the continuous rows of B as a virtual array manifold matrix from the location Then we can model the equivalent received signal as follows: where: and q represents the deterministic noise vector corresponding to the terms with the same rows as B in q.

Remarks:
The equivalent received signal exhibits a form similar to that of the true received signal expression. In spite of that, subspace-based DOA methods such as MUSIC or ESPRIT cannot be used on the equivalent received signal (say, the virtual array) directly, because p in Equation (13) behaves like a coherent source and the equivalent data are under the single snapshot condition. This, leads to a rank deficiency problem in the covariance matrix, rendering more than one source unable to be estimated and resulting in DOA estimation performance deterioration. In this paper, we tackle this problem through matrix reconstruction by using the Khatri-Rao product. Here are the definitions for subsequent derivations: e −jπ(−(N−1)r 2,y ) . . . e −jπ(−(N−1)r K,y ) e −jπ(−Nr 1,y ) e −jπ(−Nr 2,y ) . . . e −jπ(−Nr K,y ) .
We can easily get the relationship that connects B with B y and B x : D w (•), w = 1, 2, . . . , MN + M + N − 1 is a matrix operator that represents the extraction of the wth row of the matrix in the blanket and the construction of a diagonal matrix whose elements are the entries of the extracted rows. We can rewrite Equation (13) by using the following definitions: Let z w = B y D w B x p, w = 1, . . . , (MN + M + N − 1), and Equation (17) can be substituted into the following: In order to analyze the matrix reconstruction method to solve the rank problem, we ignore the noise term q in Equation (18). We notice that the Hermitian Toeplitz matrix holds the conjugate symmetric property and can be generated by a column vector whose elements locate at the non-negative positions. Inspired by that, we define: Φ y = e −jπ(N−1)r 1,y , e −jπ(N−1)r 2,y , . . . , e −jπ(N−1)r K,y ∈ C 1×K , Φ x = e −jπ(N−1)r 1,x , e −jπ(N−1)r 2,x , . . . , e −jπ(N−1)r K,x ∈ C 1×K .
After obtaining B y and B x , we can find the relations as follows: where B x (k, :) and B y (k, :) imply that the kth column of B x and B y have the dimensionality of (MN + M + N − 1) × 1, respectively. Equation (21) indicates that R B can be considered as the autocorrelation matrix of the received signal vector.
is the kth source in noiseless conditions. R B is a Hermitan and Toeplitz matrix and B x , B y are both Vandermonde matrices. After matrix reconstruction using the Khatri-Rao product, our equivalent received signal can be transferred into another virtual array received signal. The distinct elements in A B indicate that the K uncorrelated signals as defined before are incident signals on the plane array whose locations could be displaced in the set S = (x, y)|x = k x d, y = k y d 0 ≤ k x , k y ≤ MN + M + N − 2 . This newly constructed array exhibits the geometry depicted in Figure 4.

Improved Trilinear Model for Angle Estimation
Because the elevation angle   . According to the analysis in the last subsection, angle estimation can be performed as follows: It can be rewritten as a trilinear model:

Improved Trilinear Model for Angle Estimation
Because the elevation angle θ k (k = 1, 2, . . . , K) and azimuth angle φ k (k = 1, 2, . . . , K) are located at r k,x (k = 1, 2, . . . , K) and r k,y (k = 1, 2, . . . , K), we can estimate the angles in any model containing It can be rewritten as a trilinear model: Similar to Equations (17) and (18), we can obtain the following substituted expression: Equation (24) can be interpreted as slicing the three-dimensional data in a series of slices along the spatial direction [28]. The trilinear model in Equation (24) allows two more rearrangements in the matrix viewpoint, through which we obtain: x Bu = B x D u s T B y T , u = 1, 2, . . . , L, Equations (26) and (27) slice the three-dimensional data along the spatial direction and time direction, respectively. We connect the angle estimation problem of the transformed equivalent virtual array with the conventional trilinear model. The classic trilinear model-based algorithm performs the trilinear decomposition, namely PARAFAC, on the model directly. The algorithm can automatically achieve angle parameter estimation without peak searching. However, the algorithm has a heavy computational load. In order to reduce this, we propose an improved trilinear model with reduced dimensionality, which compresses the traditional model into a smaller one.
The connection between the conventional trilinear model and the improved trilinear model can be described as shown in Figure 5: The connection between the conventional trilinear model and the improved trilinear model can be described as shown in Figure 5: After compression, we can rewrite Equation (23) into: (29) After compression, we can rewrite Equation (23) into: Similarly, we easily construct the other two compressed trilinear models with respect to Equations (26) and (27):

Proposed Solution for Improved Trilinear Model
According to the three compressed trlinear models acquired from the transformation above, we can utilize the classic trilinear alternating least squares (TALS) algorithm [30] to solve the angle problem. The basic principle in TALS is using least squares (LS) to update each matrix in every step and updating other matrices in subsequent steps until algorithm convergence is achieved. The detailed steps are discussed below.
Firstly, Equation (29) gives the LS updating for matrix s : whereB y , (ŝ ) T denote for the previously obtained estimates of B y , (s ) T , respectively. y B represents the noisy compressed signal y B . Thirdly, in a similar way, Equation (31) can generate the LS solution for matrix B y :B whereB x , (ŝ ) T represents the previously acquired estimates of B x , (s ) T , respectively. x B is the noisy compressed data of x B . Finally, the matrix updating is repeated until convergence. Considering that the azimuth and elevation angles are located at B x and B y (originally in B x and B y , then compressed into B x and B y ), when we acquire the estimates of B x and B y , we need to retrieve B x and B y through Equation (28), obtaining the corresponding estimates in the form ofB x andB y .
As discussed before,B x andB y are characterized by the Vandermonde structure; we can utilize this character to solve the angle estimation problem. Assuming that (θ i , φ i ) represent the elevation and azimuth angles of the ith source, we extract the corresponding column ofB x , namely b i . We normalize the column vector b i to make the first element equal to one. After finishing those steps, the minus angle (b i ) is extracted to form a new vector: Equation (35) can be rewritten in the form of a linear equation: where Then we can obtain the LS solution: Letting D = B T B −1 B T ,we know c = Dy. Note that N is clear, so that D is a matrix of invariable dimensionality of 2 × N,which can be calculated before the signal processing. When we obtain the estimates of y, c 1 will be obtained through Equation (37) with low complexity. Similar to that, we can extract the ith column ofB y and perform the aforementioned steps to acquire the estimates of c 2 = sin φ i sin θ i . Therefore, the angle estimation results can be expressed as follows:

Complexity Analyses and Summaries
The computational complexity of our proposed method is O 4LM 2 N 2 + K + K 3 + I JFK , which is mainly caused by the covariance matrix calculation, Khatri-Rao product, and compressed trilinear decomposition.
By contrast, the method in Ref. [25] has a complexity of , where Z denotes the number of the spectral points of the total field-of-view. Due to the method of calculating two auto-correlation matrices from the decomposed two subarrays and searching for peaks in partial fields, a large computational load still exists. The method in Ref. [27] involves the step of 2D grid searching to realize the DOA estimation, which requires a complexity of O Z 3 + K 3 N 3 + M 2 L + L 3 . Therefore, our proposed method shows a superior performance in terms of computational complexity, as it does not require peak searching, or a predefined grid sampling. For the sake of clarity, the computational complexity comparison of our proposed method with other similar approaches is listed in Table 1. For spectral and grid searching with a step of 0.05 • , we also compare the complexity of methods versus snapshots and number of elements in Figures 6 and 7, respectively.

Method Complexity
Proposed The steps of the improved trilinear model-based angle estimation methods are listed in Table 2. The advantages can be summarized as follows: Firstly, we utilize the cross-correlation matrix of true received signals and extend the aperture of the CPPA in a virtual array concept, which transforms the DOA estimation problem of the CPPA into a DOA problem of the URA and augments the available DOFs in the meantime. Secondly, we reconstruct the manifold matrix of the equivalent received signal using the Khatri-Rao product and overcome the difficulty of estimating more sources under the single snapshot condition. Thirdly, we introduce the trilinear decomposition method in a compressed way into the angle estimation of the CPPA through compressing the conventional model and calculating the solutions, which reduces the computational and memory loads considerably. Finally, it's remarkable that there is no need to form predefined sampling 2D grids or search for peaks during the signal processing, which means that the algorithm can be used off-grid. the complexity of methods versus snapshots and number of elements in Figures 6 and 7, respectively.   The steps of the improved trilinear model-based angle estimation methods are listed in Table 2. The advantages can be summarized as follows: Firstly, we utilize the cross-correlation matrix of true received signals and extend the aperture of the CPPA in a virtual array concept, which respectively.   The steps of the improved trilinear model-based angle estimation methods are listed in Table 2. The advantages can be summarized as follows: Firstly, we utilize the cross-correlation matrix of true received signals and extend the aperture of the CPPA in a virtual array concept, which  Table 2. Procedure of the improved trilinear model-based angle estimation method for the CPPA.
Step 1: Obtain the equivalent received signal based on the virtual concept for the CPPA.
Step 3: Generate the conventional trilinear model of the reconstructed matrices.
Step 4: Compress the classic trilinear model with dimensionality reduction.
Step 5: Conduct the trilinear alternating least squares (TALS) algorithm with the improved model directly.
Step 6: Decompress the results in Step 5 and estimate angles via Equations (37) and (38).

Simulations and Analysis
In numerical simulations, we consider the CPPA to consist of two URAs with M = 2 and N = 3. Then the virtual array concept for the whole CPPA is a URA whose elements are located from (−2d, −2d) to (7d, 7d) in the same plane. After matrix reconstruction, the equivalent array sensors are displaced from (0, 0) to (9d, 9d). Moreover, the matrix D ∈ C 2×N in Equation (37)  in our simulation to estimate the algorithm performance. The criterion for validating the performance is the root mean square error (RMSE), and its definition is: Note that T is the number of Monte-Carlo trials, t = 1, 2, . . . , T, and K denotes the source number. φ k ,φ k denote the true azimuth angle and the estimated one, respectively. θ k ,θ k are the true elevation angle and the estimated one, respectively.
The first example is an estimation of the performance of our proposed method in terms of the signal-to-noise-ratio (SNR) and a compare of our proposed method with several DOA algorithms exploited in the CPPA, such as 2D-MUSIC and partial spectral search (PSS), which are directly performed on the array. Meanwhile, the Cramér-Rao bound (CRB) derived from Ref.
[31] for the CPPA is also employed as a benchmark to estimate the performance. Figure 8 depicts the RMSEs of different methods under the simulation conditions of T = 100 Monte-Carlo trials and an SNR in the range of −10 dB to 10 dB. The incident sources were located at (φ 1 , As shown in Figure 8, we found that the performance was improved with the increase of the SNR.
Our proposed method outperformed the other methods, exhibiting a lower complexity. Moreover, when the SNR was higher, our proposed method obtained a result closer to the CRB. The second experiment was to test the estimation performance of our proposed method versus snapshot numbers. We assumed that the SNR = 0 dB and the snapshot number varied from 50 to 500. The signals were located at          Figure 9 shows that under the condition of a small SNR and a low snapshot number, our proposed method still had a better performance than other methods and approached the CRB, owing to the aperture extension and the use of trilinear decomposition. The second experiment was to test the estimation performance of our proposed method versus snapshot numbers. We assumed that the SNR = 0 dB and the snapshot number varied from 50 to 500. The signals were located at (φ 1 , Figure 9 shows that under the condition of a small SNR and a low snapshot number, our proposed method still had a better performance than other methods and approached the CRB, owing to the aperture extension and the use of trilinear decomposition.
snapshot numbers. We assumed that the SNR = 0 dB and the snapshot number varied from 50 to 500. The signals were located at          Figure 9 shows that under the condition of a small SNR and a low snapshot number, our proposed method still had a better performance than other methods and approached the CRB, owing to the aperture extension and the use of trilinear decomposition.   It can be concluded from Figures 10 and 11 that the RMSE performance of the proposed algorithm decreases with the increase of the number of sensors in the CPPA. This is due to a larger array aperture, which can be obtained for the CPPA when more sensors are available. Moreover, no matter how many sensors were utilized, the RMSE performance continued to decrease when the SNR or snapshot number increased. This phenomenon indicates the effectiveness of our proposed  It can be concluded from Figures 10 and 11 that the RMSE performance of the proposed algorithm decreases with the increase of the number of sensors in the CPPA. This is due to a larger array aperture, which can be obtained for the CPPA when more sensors are available. Moreover, no matter how many sensors were utilized, the RMSE performance continued to decrease when the SNR or snapshot number increased. This phenomenon indicates the effectiveness of our proposed method for different CPPA configurations. When we consider real applications, such as the It can be concluded from Figures 10 and 11 that the RMSE performance of the proposed algorithm decreases with the increase of the number of sensors in the CPPA. This is due to a larger array aperture, which can be obtained for the CPPA when more sensors are available. Moreover, no matter how many sensors were utilized, the RMSE performance continued to decrease when the SNR or snapshot number increased. This phenomenon indicates the effectiveness of our proposed method for different CPPA configurations. When we consider real applications, such as the measurement of the radar signal parameters, we only need a small sensor number to construct the co-prime planar array under a large snapshot number condition. Figure 11 shows that the array figuration of M = 2, N = 3 can achieve a similar angle estimation accuracy to the figuration of M = 4, N = 5 when the snapshot number is 500. This result means that we can utilize limited sensors to obtain high accuracy, which is also beneficial for minimizing array arrangement costs.
Through the above three experiments, we validated the effectiveness of our method used in co-prime planar arrays and compared it with similar state-of-the-art approaches from the perspective of accuracy and sensor numbers. It is known that high accuracy and low figuration sensor array cost are very important in complex electromagnetic environments where multiple effective signals are submerged by various types of noises and interferences. Even in low SNR environments (such as SNR = −15 dB), the proposed method can still detect the target signals and reduce the side lobes of the beam patterns. Figures 12 and 13 show the results of the angle estimation at the desired angle (azimuth = 30 • , elevation = 45 • ). As shown in Figures 12 and 13, the side lobes are lower than the direct synthesis beam patterns and the main lobe is nearly identical to them.

Conclusions
In this paper, a virtual array concept is considered by utilizing the cross-correlation matrix of the received data for a CPPA, which increases the available DOFs to a relative degree. After

Conclusions
In this paper, a virtual array concept is considered by utilizing the cross-correlation matrix of the received data for a CPPA, which increases the available DOFs to a relative degree. After obtaining the virtual array, matrix reconstruction using the Khatri-Rao product was proposed to

Conclusions
In this paper, a virtual array concept is considered by utilizing the cross-correlation matrix of the received data for a CPPA, which increases the available DOFs to a relative degree. After obtaining the virtual array, matrix reconstruction using the Khatri-Rao product was proposed to tackle the problem of matrix rank deficiency and to ensure that the trilinear decomposition method can be introduced to estimate more sources under finite snapshots. Furthermore, compressed matrices were applied before the trilinear decomposition, which constructs an improved trilinear model to be solved. This step dramatically reduces the computational complexity and memory loads for the angle estimation, rendering the algorithm easier to implement, to some extent. More importantly, our proposed method requires no peak searching nor grid sampling, and achieves automatic pairing. This means that the proposed method can achieve high estimation accuracy with low complexity compared to other similar approaches. Furthermore, our proposed method only needs finite sensor numbers to implement the source number detection tasks, which efficiently reduces the cost of sensor array arrangement in real conditions.