Two-Dimensional Direction-of-Arrival and Polarization Parameter Estimation Using Parallel Co-Prime Polarization Sensitive Array

Target detection is critical in many mission critical sensors and sensor network (MC-SSN) applications. For target detection in complicated electromagnetic environment, DOA estimation using polarization sensitive array (PSA) has been receiving increased attentions. In this paper, we propose the parallel co-prime polarization sensitive array (PCP-PSA) which consists of the cocentered orthogonal dipole triads (CODTs) to estimate two-dimensional direction-of-arrival (2D DOA) and polarization parameters. The degrees of freedom (DOFs) have been extended due to the co-prime structure, so that the more signals can be detected and the estimation accuracy is improved. In order to reduce the computation complexity, we construct a new cross-covariance matrix based on the CODTs, which converts the two-dimensional DOA estimation into two independent one-dimensional DOA estimations. Then, the spatial smoothing-based multiple signal classification algorithm(MUSIC) and the sparse representation-based method are applied to estimate 2D DOA with only one-dimensional (1D) peak searching and 1D dictionary, respectively. Finally, the polarization parameters are estimated by using the cross-covariance matrix between components of electric field vector. Compared with previous PSA-based algorithms, the proposed algorithm based on PCP-PSA can solve the underdetermined 2D DOA and polarization parameters estimation problem and has better estimation accuracy. Theoretical analyses and simulation results verify the effectiveness of the proposed methods in terms of computational complexity and estimation accuracy.


I. INTRODUCTION
In recent years, mission critical sensors and sensor network (MC-SSN) applications, such as target detection and reliable communication, have received increasing attentions from both research community and industry [1]- [3]. The sensors network composed of the scalar sensors used to target detection on radar with DOA information. Moreover, in industry [4], the series of the sensors used to construct the sensors array and networks to get DOA information and establish communication [5], [6].DOA information [7] is critical in these MC-SSN applications. Compared with the scalar The associate editor coordinating the review of this manuscript and approving it for publication was Qilian Liang . sensor array, DOA estimation using the polarization sensitive array (PSA) offers better performance in estimation accuracy, recognition accuracy and anti-jamming [8]- [10]. Therefore, various effective DOA estimation algorithms based on PSA have been proposed [11]- [13]. The Multiple Signal Classification algorithms(MUSIC) and Estimation of Signal Parameters via Rotational Invariance Techniques(ESPRIT) were transferred to PSA. The polarized MUSIC with multi-dimensional peak searching of DOA and polarization parameters have been proposed with a long-vector (LV) data model of the electromagnetic wave signal [13], [14]. A series of polarized ESPRIT algorithms [15]- [17] based on time, space, polarization invariance were presented to reduce the computational complexity. The vector cross-product based algorithm [18], [19] can estimate the parameters without ambiguity, where only a single six-component vector-sensor is used. In general, algorithm using PSA can achieve better performance with the vector structure of electromagnetic signal.
As the modern electromagnetic environment is becoming more complex, the increasing number of the radiation sources and wide distribution density is become a common situation. In this situation, the number of incident signals usually is unknown, and underdetermined (the number of signals is more than the number of array elements) DOA estimation problems are becoming more and more common. Therefore, research on under-determined DOA estimation has important practical significance.
The aforementioned algorithms can not solve the under-determined problem that the number of signals is more than array elements, and the inter-element spacing limit of less than or equal to half a wavelength also restricts the estimation performance [20], [21]. In order to solve these problems, coprime array is proposed as a new sparse array technology, which provides a new idea for solving the underdetermined DOA estimation problem. Co-prime array (CPA) [22]- [25] was presented to improve the degrees of freedom (DOFs) and estimation accuracy. By vectorising auto-covariance matrix array of the received data, the DOFs of the CPA could reach O(MN ) with only O(M + N ) physical sensors when m and n are co-prime numbers. M and N separately express the number of two subarrays' array elements in the CPA. The parallel co-prime array (PCPA) [25]- [27], three-parallel co-prime array (TPCPA) [28], and co-prime planar array (CPPA) [20]- [21], [29], [30] are constructed by vectorising the cross-covariance matrix of multiple subarrays for 2D DOA estimation. In order to solve the correlation during the vectorisation, the spatial smoothing(SS) technique based MUSIC [31]- [37] was introduced to the CPA. Subsequently, the sparse representation (SR) framework based CPA algorithms [37]- [42] were developed to extend the DOFs. In summary, setting up the sensors array with the one or two dimensional co-prime structure can improve the DOFs, thereby more signals can be detected. Meanwhile, the array aperture is extended because of the sparse configuration of CPA, which improves the estimation performance. However, the CPA has not been introduced to the polarization sensitive array models.
In this paper, we propose a parallel co-prime polarization sensitive array (PCP-PSA), and construct a novel cross-covariance matrix. Compared with the existing estimation algorithm based on PSA, the proposed model dramatically increases the DOFs since the co-prime structure is constructed, which is applicable to underdetermined situations and improves the estimation performance. More importantly, in order to reduce the computational complexity, the four parameters estimation problem (2D DOA and 2D polarization parameters) is decoupled into two 1D-DOA parameter and one 2D parameters parameter estimation. Through the vectorisation operation, a virtual uniform linear array (ULA) with M (N + 1) DOFs is constructed. Then proposing separately using SS based MUSIC or SR based method in sequence to estimate the angle of 2D DOA. The simulation results indicate that the proposed algorithm has superior performance on parameter estimation than existing PSA based estimation algorithm.
The remainder of this paper is organized as follows. The problem is formulated in Section 2. The DOA estimation and polarization parameter estimation of the proposed algorithm are explicitly described in Section 3, respectively. The numerical simulations are conducted to validate the effectiveness of the proposed algorithm in Section 4. Finally, the paper is concluded in Section 5.

II. PROBLEM FORMULATION A. PARALLEL CO-PRIME POLARIZATION SENSITIVE ARRAY
Consider a parallel co-prime polarization sensitive array that consists of two uniform linear subarrays in the xoy plane, as shown in Figure 1. The elements here in are all the cocentered orthogonal dipole triads (CODTs) whose three dipoles parallel to x−, y− and z−axis, respectively.The subarray 1 has 2M sensors with equal interelement spacing Nd 1 ,whereas the subarray 2 has N sensors with equal interelement spacing Md 1 , where d 1 = λ/2 is a fundamental spacing, and λ is the signal wavelength. The displacement spacing between the two subarrays is d 2 = λ/2. The numbers M and N are co-prime, without loss of generality, assume that M < N . Obviously, the PCP-PSA contains 2M + N physical sensors.

B. 2D DOA ANGLES DEFINITION
In the three-dimensional space, a pair of two-dimensional angles is required to characterize the incident direction of the electromagnetic wave signal. As illustrated in Figure 2, θ k ∈ [0, 2π) and φ k ∈ [0, π/2] denote the azimuth angle measured from the positive x−axis to the projection of the incident direction on the xoy plane and the elevation angle measured from the vertical z−axis, respectively, whereas α ∈ [0, π] and β ∈ [0, π] are spatial angles measured from the positive x−axis and the positive y−axis, respectively. These two pairs of two-dimensional spatial angles both can be used to determine the DOA of incident signal, and have the following relationship.
1) The received signal data is statistically independent among the sensors, dipoles and snapshots.
2) The complex noise is supposed to be temporally and spatially white Gaussian and uncorrelated with the sources.
3) There is no mutual coupling effect among the sensors and dipoles.

III. PARAMETER ESTIMATION A. TWO-DIMENSIONAL DOA ESTIMATION
The second-order statistic of the array received data is utilized to estimate the parameters. The cross-covariance matrix of x-componets of two subarrays is where 2 ) is the signal cross-covariance matrix of x−componets, and σ 2 k denotes the power of k−th signal. The superscript (·) H denotes the complex conjugate transpose operation, E[·]denotes the expectation operation. Similarly, the cross-covariance matrices of y−componet and z−componet can be obtained as where S y = diag(σ 2 1 e y,1 2 , σ 2 2 e y,2 2 , · · · , σ 2 K e y,K 2 ) and . It can be seen that there are no noise terms in the cross-covariance matrices, hence it would less affected by noise. According to (2), the following quality can be obtained Therefore, a new cross-covariance matrix can be constructed as where S = diag(σ 2 1 , · · · , σ 2 K ) is the covariance matrix of signals. It is obviously that R only involves DOA angles but not the polarization parameters. According to the relationship between (θ, φ) and (α, β), the new cross-covariance matrix R can also be written as R = A [1] (β)S A [2] where a [1] (β k ) = [1, e j2πNd 1 cos β k /λ , · · · , e j2π(2M −1)Nd 1 cos β k /λ ] T and a [2] (β k ) = [1, e j2πMd 1 cos β k /λ , · · · , e j2π(N −1)Md 1 cos β k /λ ] T denote the k-th column of A [1] (β) and A [2] (β), respectively.
It is worthy to note that the 2D DOA angles (α, β) have been separated into two independent parts in (10), and the angle α has been transformed into the diagonal entries of the matrixS.
By using the vectorisation operator, the new crosscovariance matrix R can be transformed into a vector form as (11) where and (·) * denote the Khatri-Rao product and conjugate operation, respective. s is a column vector that consists of the diagonal entries ofS, b(β k ) = [a [2] (β k ) ⊗ a [1] where ⊗ denotes the Kronecker product. The vector r can be regard as the received signal data from a virtual linear array with the sensor locations at where x [1] i is the location at x−axis of the i−th element of subarray 1, and x [2] j is the location at x−axis of the j-th element of subarray 2. s is the single-snapshot signal data, and B(β) is the corresponding steering matrix. Benefiting from the co-prime characteristic of numbers M and N , we can construct a virtual ULA with sensor locations at By extracting the entries and rows corresponding to the locations inL from r and B(β), the single-snapshot received data of virtual ULA can be given by . It can be seen that the DOFs of the PCP-PSA are extended from 2M + N to M (N + 1), which allows that the DOA angles can be estimated under the underdetermined condition.
Since the vectorisation operation in (11) induces correlation between the incident signals, hence, the signal model in (14) would be considered as K 'correlated' incident signals impinging on the virtual ULA. Due to the correlation of signals and only single available snapshot in (14), the subspacebased DOA estimation algorithm cannot be applied directly.
The spatial smoothing method can effectively recover the rank of the array covariance matrix, so that the correlated signals can be estimated by subspace-based algorithm. On the other hand, the sparse representation based DOA estimation algorithms can handle the correlated problem naturally, and can work properly even if the snapshot is insufficient. In the following, the forward/backward spatial smoothing (FBSS)-MUSIC algorithm [32] and 1−norm penalty-based sparse recover algorithm are used to estimate the DOA angle β.
According to (14), the virtual ULA with M (N + 1) sensors can be divided into P 1 overlapping subarrays with P 2 sensors in each subarray, and P 1 ,P 2 satisfy P 1 + P 2 = M (N + 1) + 1. If the subarrays are arranged with forward spatial smoothing, the first subarray consists of the sensors located at 0 to (P 2 − 1)d 1 . The received signal data of the p−th subarray is where p = 1, 2, · · · , P 1 . Then, the forward spatial smoothing covariance matrix of the virtual ULA can be calculated by Similarly, if the subarrays are arranged with backward spatial smoothing, the first subarray consists of the sensors located at P 1 d 1 to (P 1 + P 2 − 1)d 1 . The received signal data of the p−th subarray is where p = 1, 2, · · · , P 1 . Then, the backward spatial smoothing covariance matrix of the virtual uniform linear array can be calculated by where  is a P 2 × P 2 permute matrix with back-diagonal of 1.
Thus the FBSS covariance matrix can be obtained as Appling the FBSS covariance matrix to the conventional MUSIC algorithm, the DOA angle β can be estimated. Note that the FBSS covariance matrix is a full-rank P 2 ×P 2 matrix, hence, it is can be used to estimate the DOA angle β for P 2 − 1 incident signals. According to the theory of spatial smoothing, the element number in one subarray and the subarray number are mutually restricted, and estimation performance would degrade when either of them is too large or too small. Then the sparse representation based DOA estimation algorithm is discussed. By discretizing the whole β angle domain, a sampling gridβ 1 ,β 2 , · · · ,β Q with K is formed. Assume the grid is dense enough so that the actual DOA angles β k ( k = 1) K only lie within the grids. To obtain denser grids with less computational complexity, the multiresolution grid refinement [43] can be used. Then the data vector z in (14) can be sparsely represented as VOLUME 8, 2020 where φ(β) = [b(β 1 ),b(β 2 ), . . . ,b(β Q )] ∈ C (M (N +1)×Q) is an overcomplete dictionary, and u = [u 1 , u 2 , · · · , u Q ] T is a K −sparse coefficient vector with u q = 0 if ∃β q = β k and u q = 0 otherwise. It can be seen that φ(β) only depends on the DOA angle β, whereas u depends on the other DOA angle α and the signal power. Hence, the estimation of DOA angle β can be obtained by solving the following optimization problem where · 2 and · 1 denote the 2 − and 1 −norm, respectively. µ is the regularization parameter that balances the trade-off between the residual error of z and the sparsity of u. In fact, the above optimization problem is a second-order cone program (SOCP) problem [44], and can be efficiently solved by off-the-shelf optimization software package such as CVX. Therefore, the DOA estimation problem turns out to be that of recovering theK sparse vector u and detecting the locations of nonzero elements therein. Once the estimation of DOA angle β is obtained, the steering matrixB(β) of the virtual ULA in (14) can be calculated aŝ where {β k } K k=1 is the estimation results. Then, the estimation of vector s can be calculated by the least squares method asŝ where (·) −1 represents the matrix inverse operation. Therefore, the estimation of {α K } K k=1 can be obtained bŷ α k = cos −1 (−λ arg(ŝ k )/2πd 2 ), whereŝ k is the k-th entry ofŝ, and arg(·) represents the argument. It is obviously that the DOA angles {α k } K k=1 and {β k } K k=1 are paired automatically.

B. POLARIZATION PARAMETER ESTIMATION
In this section, we estimate the polarization parameters based on the DOA estimation results. From (5) to (7), the crosscovariance matrices between the x−, y−components of subarray 1 and z− component of subarray 2 can be constructed as where S xz = diag(σ 2 1 e x,1 e * z,1 , · · · , σ 2 k e x,K e * z,K ) and S yz = diag(σ 2 1 e y,1 e * z,1 , · · · , σ 2 K e y,K e * z,K ) both are diagonal matrices. By using the vectorisation operator, we can get [2] (θ, φ) * A [1] (θ, φ) · s yz , (28) where s xz ands yz are column vectors consist of the diagonal elements of S xz and S yz , respectively. Since 2D DOA angle estimation results are obtained, s x z and s y z can be estimated by LS method asŝ whereÃ = A [2] (θ ,φ) * A [1] (θ ,φ) can be calculated with {θ k ,φ k } K (k=1) . Therefore, the estimation of polarization parameters can be obtained bŷ whereŝ xz,k andŝ yz,k are the k−th entries ofŝ xz andŝ yz , respectively. σ 2 k = ŝ k 2 is the power of k-th signal which can be obtained from the results of (23). It can be seen that the polarization parameter estimates {γ k ,η k } K k=1 and DOA angle estimates {θ k ,φ k } K k=1 are also paired automatically in the process of the least square operation.

IV. SIMULATION RESULTS
In this section, a series of numerical simulations under different conditions are conducted to investigate the estimation performance of the proposed algorithm. The results are compared with the PCPA algorithm [20], TPCPA algorithm [23], and LV-MUSIC algorithm [7]. The parallel co-prime array is used in the proposed algorithm and PCPA algorithm, and the co-prime numbers are M 1 = 5 and N 1 = 7, hence, there are 17 physical sensors and 40 DOFs. The there-parallel coprime array is used in the TPCPA algorithm, and co-prime numbers are M 1 = 4 and N 2 = 5, hence, there are 18 physical sensors and 47 DOFs. The array used in the LV-MUSIC is same as the proposed algorithm. Moreover, the array elements are CODTs for the proposed algorithm and LV-MUSIC algorithm, whereas they are scalar sensors for PCPA and TPCPA algorithm.

A. THE SPATIAL SPECTRUMS UNDER THE UNDERDETERMINED CONDITION
In this subsection, we compare the spatial spectrums of the DOA angle β under the underdetermined condition. Assume there are K = 18 far-field uncorrelated narrowband signals impinging on the array. The DOA angle β are distributed within 18 • to 154 • with step of 8 • , the DOA angle α are distributed within 13 • to 77 • with step of 8 • , the auxiliary polarization angle γ and the polarization phase difference η are randomly distributed within the ranges of (−90 • ,90 • ) and (−180 • ,180 • ), respectively. The snapshot number is 500 and the SNR is 10dB.
We firstly compare the spatial spectrum of the proposed algorithm based on SR with that of the PCPA algorithm. The uniform spatial grid in β angle domain is formed with interval of 0.1 • . The regularization parameter µ is 0.3 for the proposed algorithm and 1 for PCPA algorithm. As it can be seen in Figure 3, both these two algorithms can obtained the correct estimation of β under the underdetermined condition. However, the PCPA algorithm has larger estimated bias than the proposed algorithm, and there are visible pseudo-peaks in the spatial spectrum. This certifies that the proposed algorithm has better estimation stability than PCPA algorithm. Then, we compare the spatial spectrums of the proposed algorithm based on FBSS-MUSIC with different sensor number of each overlapping subarray, as shown in Figure 4. It is obviously that sensor number of subarrays seriously affects the spatial spectrum. When the subarrays contain fewer sensors, the deviations of some peaks increase since more entries in the original covariance matrix R are missed. For example, there are some accuracy loss when the DOA angle β are 18 • , 26 • , 146 • , 154 • . This is mainly because the spatial smoothing only uses the diagonal area elements of the covariance matrix of the virtual uniform linear array. In this situation, when the subarrays contain fewer sensors, more non diagonal area elements will be lost. When the subarray contains more sensors, the orthogonality of signal-and noise-subspace would decrease since the spatial smoothing time reduces. Therefore, as previously discussed, the moderate choice of subarray sensor number is significant.

B. ESTIMATION PERFORMANCE VERSUS SNR AND SNAPSHOT NUMBER
In this section, we discuss the 2D DOA and polarization estimation accuracy versus SNR and snapshot number. Assume there are K = 4 far-field narrowband signals whose DOA angels and polarization parameters The entire spatial domain is considered, and the searching step or grid interval is 0.1 •• . When the snapshot number is fixed to 500, the SNR varies from 5dB to 30dB with the step size 5dB; when the SNR is fixed to 10dB, the snapshot number varies from 100 to 1000 with the step size 100. Figure 5 and 6 show the root mean squared errors (RMSEs) of 2D DOA estiamtion versus SNR and snapshot number. It can be seen that the RMSEs of DOA for these algorithms all decrease with the increase of SNR and snapshots. Due to the reconstructed cross-covariance matrix, the proposed algorithm outperforms than the other algorithms. In the case of low SNR and few snapshots, the orthogonality of the signal-and noise-subspace in the MUSIC-based algorithms is affected, leading to a slight increase of the estimation RMSE than SR-based algorithms. Moreover, since the entries in the covariance matrix of data vector z would be missed as we discussed in the first simulation, the proposed algorithm based on SR has more accurate estimate results than that based on FBSS-MUSIC.
Next, Figures 7 and 8 show the estimation accuracy of the polarization parameters versus SNR and snapshot number. Similar to the RMSE curves of DOA estimation, the RMSEs of polarization parameters decreases with the increase of SNR and snapshots. As shown in Figure 7, in the case of low SNR, the proposed algorithm has more accurate estimation results since the usage of cross-covariance matrix which is less affected by noise. When the SNR is high, LV-MUSIC VOLUME 8, 2020  algorithm have smaller estimate RMSEs of polarization parameters than the proposed algorithm based on FBSS-MUSIC. The main reason is that the proposed algorithm estimates the polarization parameters based on DOA estimation results. The DOA estimation deviation would cause more estimation bias in the process of polarization parameter estimation. This would not appear in the LV-MUSIC algorithm. However, the proposed algorithm based on SR has enough accurate DOA estimation, so that it has the smallest polarization parameter estimation RMSE.   500 snapshots. We compare the proposed algorithm with the LV-MUSIC algorithm and PCPA algorithm, then their algorithm complexity analysis and average running time are shown in Table 1. It can be seen that the running time of the proposed algorithm based on FBSS-MUSIC is drastically reduced compared with LV-MUSIC algorithm. This is mainly because the proposed algorithm constructs the virtual ULA, and only 1D peak searching is required to estimate 2D DOA. On the other hand, since the spectrum function of the proposed algorithm is simpler than LV-MUSIC algorithm, the computational amount of spectrum function on each searching point is also reduced. In addition, the proposed algorithm based on SR has the similar running time with PCPA algorithm, which is consistent with the aforementioned analysis of theoretical computational complexity.
In this subsection, we discuss the major computational complexity of the proposed algorithm. The construction of virtual ULA and its data vector z requires O(MNL) flops. In the algorithm based on FBSS, comparing with the Kroenke product in the algorithm, the rest of algorithm complexity is negligible. Therefore, the complexity of the algorithm mainly comes from the operation of the matrix. In addition, in the algorithm based on the SR. When the number of array elements is the same, the algorithm complexity of parameter estimation algorithm under the sparse representation is huge, and it is usually several times more expensive than DOA parameters estimations. Therefore, the dictionary mainly determines the algorithm complexity of the algorithm based on the SR. The estimation of spatial angle β based on FBSS-MUSIC or sparse recover requires O(N β P 3 2 ) or O(Q 3 ) flops, respectively, where N β denotes the searching point number in angle β domain; the estimation of spatial angle α requires O(M (N + 1)K 2 ). Typically assuming N β , Q, L P 2 > K , M , N , therefore the computational complexity of the proposed algorithm is mainly in the 1D peak searching process of FBSS-MUSIC or solving the optimization problem in (21). By comparison, the current 2D MUSIC-like algorithms need 2D peak searching process, such as LV-MUSIC algorithm [7], which requires O(N θ N φ (2M + N ) 3 ) flops where N θ and N φ denote the searching point number in azimuth and elevation domain. This would lead much more computational amount than the proposed algorithm based on FBSS-MUSIC. On the other hand, 2D DOA estimation based sparse recover all need to solve the optimization problem which requires the similar computational complexity as the proposed algorithm based on SR.

V. CONCLUSION
In this paper, a novel 2D DOA and polarization parameter estimation algorithm based on parallel co-prime polarization sensitive array is proposed. We construct a new cross-covariance matrix that does not contain polarization parameters, thus the 2D DOA and polarization parameters can be estimated separately. Through a vectorisation operation, the virtual ULA is constructed. Therefore, the DOFs of the array increase and more signals can be detected with the limited number of physical sensors. Then, the FBSS-MUSIC with 1D peak searching with 1D dictionary and SR-based method are used for DOA estimation. The FBSS-MUSIC method needs to construct overlapping subarray to solve the correlation of signals, thus causing the loss of DOFs. In contrast, the SR based method can naturally process the correlated signals and fully utilize the DOFs, achieving the detection of more signals and higher estimation accuracy. Since the co-prime structure is applied to PSA, the DOFs of the array are extended. Meanwhile, due to the sparse placement of sensors, the array aperture is increased and the mutual coupling is weakened. From the simulation results, it can be seen that these advantages make the proposed algorithm can solve the underdetermined estimation problem and compared with PCPA algorithm, PCP-PSA algorithm provides better estimation stability in underdetermined conditions.