A PARAFAC-based algorithm for multidimensional parameter estimation in polarimetric bistatic MIMO radar

In this article, we investigate the problem of applying the parallel factor quadrilinear decomposition technique to multidimensional target parameter estimation in a polarimetric bistatic multiple-input multiple-output (MIMO) radar system with a uniform rectangular array at the transmitter and a cross-dipole-based uniform rectangular array at the receiver. The signal model is developed, and a novel algorithm is proposed exploiting the quadrilinear alternating least squares to jointly estimate the two-dimensional direction of departure (2D-DOD), two-dimensional direction of arrival (2D-DOA), polarization parameters and Doppler frequency. Multidimensional parameters can be automatically paired by this algorithm to avoid the performance degradation resulting from wrong pairing. The developed algorithm requires neither multidimensional spectral peak searching nor covariance matrix estimation and several eigen-value decompositions that may bring error accumulation. Furthermore, multiple targets having close 2D-DODs and close 2D-DOAs or even the same 2D-DOD or 2D-DOA are distinguishable by means of polarization diversity. The algorithm improves the performance of multi-target identification and three-dimensional localization. Numerical simulations demonstrate the effectiveness of the proposed algorithm.


Introduction
In recent years, joint parameter estimation in multipleinput multiple-output (MIMO) radar [1][2][3] has drawn considerable attention for target identification, localization, imaging, etc. Specifically, bistatic MIMO radar with respectively colocated transmitter and receiver array antennas is widely investigated by many researchers for its capability of jointly estimating the direction of departure (DOD) and direction of arrival (DOA) of targets. Many spatial spectrum estimation approaches have been developed for joint DOD and DOA estimation in a bistatic MIMO radar [4][5][6][7][8][9][10][11]. In [4], a two-dimensional (2D) Capon-based method was developed by searching through all the 2D space to find the DOA and DOD of a target. By exploiting the invariance property of both the transmitter and receiver arrays, some *Correspondence: jiangh@jlu.edu.cn College of Communication Engineering, Jilin University, Changchun 130012, China approaches using the estimation of signal parameter via rotational invariance techniques (ESPRIT) [5][6][7] have been presented for joint DOA and DOD estimation, avoiding the exhaustive peak searching. However, an additional pair matching between the DOAs and DODs is required for ESPRIT-like methods. Besides, polynomial root finding [8] and combined ESPRIT-multiple signal classification (MUSIC) approach [9] were investigated by decomposing the 2D angle estimation into double one-dimensional (1D) angle estimation, allowing an automatic pairing between the DOAs and DODs. To reduce computation load, a reduced-dimension MUSIC method was proposed in [10]. Considering the characteristics of non-stationary target signals, joint estimation of DOD and DOA information of maneuvering targets was examined in [11] exploiting spatial timefrequency distribution. For moving targets, the literatures in [12][13][14] addressed to joint DOD, DOA and Doppler frequency estimation. http://asp.eurasipjournals.com/content/2013 /1/133 Notably, the polarization state of a target will change upon reflection. Polarization diversity has been proved to be important in target identification especially when multiple targets are so closely spaced that they cannot be distinguished well in spatial domain. In [15,16], joint DOD-DOA-polarization estimation and joint DOD, 2D-DOA and polarization estimation exploiting the ESPRIT technique were proposed for bistatic MIMO radar. The simulations have shown that multiple targets having close DODs or DOAs but different polarizations own different array manifolds and are distinguishable based on their polarization diversity.
Recently, parallel factor (PARAFAC) analysis, as an analysis method of high-dimensional data, has become a new technology applied to signal processing and communication field [17,18]. The parallel factor model is a generalization of low-rank decomposition to three-or multi-way arrays, which was first introduced as a data analysis tool in psychometrics. Over the past decades, PARAFAC ideas have been applied in multiple-invariance sensor array processing with emphasis on identifiability results [18]. In recent years, PARAFAC has become a new research means in MIMO radar [19][20][21]. The PARAFAC analysis algorithms [19,20] and adaptive PARAFAC algorithm [21] have been developed for the estimation of DOAs and DODs of multiple targets.
Despite the fact that several PARAFAC-based methods have been proposed for DOA and DOD estimation in bistatic MIMO radars, trilinear decomposition algorithms [17][18][19][20][21] with trilinear alternating least square (TALS) are commonly used to fit PARAFAC model. In this article, we investigate joint estimation of DOD, DOA, Doppler frequency and polarization. The polarization can be considered as the fourth dimension and hence a four-way description happens to fit better as long as the targets have diverse polarizations. Thus, PARAFAC quadrilinear decomposition [22,23] with quadrilinear alternating least square (QALS) is applied to multidimensional parameter estimation in polarimetric bistatic MIMO radar.
With the restriction that targets can only be located on a 2D plane in previous algorithms based on 1D-DOD and 1D-DOA estimation using uniform linear arrays, the estimation of 2D-DOD and 2D-DOA is investigated in this article exploiting a uniform rectangular array at the transmitter and a cross-dipole-based uniform rectangular array at the receiver. The signal model for polarimetric bistatic MIMO radar is developed, and a novel PARAFAC QALS-based algorithm is presented for joint estimation of seven target parameters, including the 2D-DOA, 2D-DOD, two polarization parameters and Doppler frequency. The proposed algorithm avoids multidimensional spectral peak searching, covariance matrix estimation and several eigen-decompositions that may bring error accumulation, which can enhances the accuracy of estimation. Multidimensional parameter pairing is obtained automatically by this algorithm, which can avoid the performance degradation resulting from wrong pairing. Furthermore, polarization diversity of multiple target characteristics is exploited to distinguish multiple targets having close 2D-DODs and 2D-DOAs or even the same 2D-DOD or 2D-DOA, which can improve the resolution of multi-target identification and 3D localization. The merits of the algorithm are investigated via numerical simulations.
This paper is organized as follows: in the next section, the signal model for polarimetric bistatic MIMO radar is developed. A novel PARAFAC quadrilinear decomposition-based algorithm for 2D-DOA, 2D-DOD, polarizations and Doppler frequency estimation is presented in section 3. In section 4, the results of simulation are given to verify the performance of the proposed method. Finally, a conclusion is drawn in section 5.

Received signal
As illustrated in Figure 1, consider a polarimetric bistatic MIMO radar equipped with a uniform rectangular transmit array having M = M 1 M 2 sensors and a uniform rectangular receive array having N = N 1 N 2 crossed dipole sensors, where M 1 ,M 2 and N 1 ,N 2 are the numbers of the transmit and receive sensors in the x-axis and y-axis, respectively. At the receiver array, the crossed dipoles are used to measure the polarization states of a transverse electromagnetic (TEM) wave and exploit polarization diversity. The antennas are of ideal, identical isotropic sensors. The inter-sensor spacings at the transmit and receive arrays are d t and d r , respectively, their lengths being no more than half wavelength (λ/2). At the transmit site, M different coded signals with identical bandwidth and center frequency are emitted simultaneously, which are temporally orthogonal. Denote s m = [ s m (1), · · · , s m (K)] the mth transmit waveform vector having a code length of K and a unit norm, m = 1, · · · , M. There are P targets located in the same range-bin of the receiver site. They are assumed to be point sources. For the pth target, p = 1, · · · , P, we denote θ tp and φ tp its transmit elevation and azimuth angles (namely, 2D-DOD), θ rp and φ rp its receive elevation and azimuth angles (namely, 2D-DOA), γ p and η p its polarization angle and polarization phase difference, and f dp its Doppler frequency.
Denote X (l) m,n,p the reflected signal of the pth target from the mth transmit sensor to the nth receive sensor in the lth snapshot, l = 1, · · · , L, m = 1, · · · , M, n = 1, · · · , N, p = 1, · · · , P, and X (l) 1,1,p the reflected signal of the pth target from the reference transmit sensor to the reference http://asp.eurasipjournals.com/content/2013/1/133 target y ' t t r Transmit 1 2 Receive array receive sensor. The incoming signals of fully polarized electric field TEM wave impinge on the reference receive sensor with crossed dipole antenna. For the pth target, we assume that it arrives at the receive array with the elevation/azimuth angles θ rp and φ rp . Thus, the observed signal X (l) 1,1,p consists of two polarization components in the x direction and y direction [24], the composite of which can be written as 1,1,p = cos θ rp cos φ rp − sin φ rp cos θ rp sin φ rp cos φ rp sin γ p e jη p cos γ p β p e j2π f dp t l where β p denotes the complex amplitude proportional to the radar cross section of the pth target. The target amplitude remains constant in the L snapshots. t l = lT s is the slow time with T s being the duration of a snapshot. Z (l) 1,1,p is the noise measured at the reference sensor. Since the observed signal X (l) m,n,p coming from the mth transmit sensor to the nth receive sensor has the phase shifts relative to X m = 1, · · · , M, n = 1, · · · , N, p = 1, · · · , P, m i = 1, · · · , M i , n j = 1, · · · , N j , i, j = 1, 2, and m,n,p is the observed noise. Collect the data of X (l) m,n,p for m = 1, · · · , M, n = 1, · · · , N, p = 1, · · · , P, forming a 2N × K complex matrix where ⊗ denotes the Kronecker product, and is a 2 × 1 polarization vector for a crossed dipole antenna array, 0 ≤ γ p ≤ π/2, −π ≤ η p < π. The pth steering vectors with respect to the transmit and receive uniform rectangular arrays are written as a t (θ tp , φ tp ) ∈ C M×1 and a r (θ rp , φ rp ) ∈ C N×1 , respectively, with where The observed matrix in (4) can be rewritten as where represents the Khatri-Rao product (columnwise Kronecker product). A t (θ t , φ t ), A r (θ r , φ r ) and G(θ r , φ r , γ , η) denote the transmit steering matrix, receive steering matrix and polarization matrix,

Pulse compression and vectorization
The received signals are pulse-compressed using M transmit waveforms, forming the output matrix X (l) = X (l) S H , which can be expressed by where Z (l) S H denotes the noise matrix after pulse compression. Vectorizing X (l) in (9) yields the following 2MN column vector: where A(θ t , φ t , θ r , φ r , γ , η) denotes the total steering matrix with In (10) and (11), the transformation of vec(ABC) = (C T ⊗ A)vec(B) has been used. Collecting all the observed data of y (l) over L time slots yields the following 2MN × L complex matrix with Y = y (1) , · · · , y (L) and N = n (1) , · · · , n (L) representing the observed matrix and noise matrix. B = b (1) , · · · , b (L) T denotes the target matrix, which contains the estimated parameter of Doppler frequency f d .
Based on the signal model of polarimetric MIMO radar developed in (12), the problem of interest is to jointly estimate the multidimensional parameters θ t , φ t , θ r , φ r , γ , η and f d .

PARAFAC quadrilinear decomposition-based algorithm for multidimensional parameter estimation
In this section, we present the multidimensional parameter estimation algorithm in polarimetric MIMO radar using PARAFAC quadrilinear decomposition.

Definition 1. (Quadrilinear decomposition in tensor format) A quadrilinear decomposition of a four-order ten
and X (4) ∈ C JKL×I be the four matrix representations of a four-way array. The quadrilinear decomposition of X can be written as four equivalent matrices X (1) Under certain conditions, X can be decomposed uniquely into matrices A, B, C and D. These conditions are based on the notion of Kruskal-rank [17,18]. Definition 3. (Kruskal-rank or k-rank) [18] Consider a matrix A ∈ C I×M . If rank(A) = r, then A contains a collection of r linearly independent columns. Moreover, if every k < M columns of A are linearly independent, but this does not hold for every k + 1 columns, then A has k-rank k A = k. Note that k A < rank(A), ∀A. Theorem 1. (Uniqueness of quadrilinear decomposition) [25] Consider that a matrix representation of X is X = (A B C) D T , where A ∈ C I×M , B ∈ C J×M , C ∈ http://asp.eurasipjournals.com/content/2013/1/133 C K×M and D ∈ C L×M are the four mode matrices of X, and M denotes the common dimension. If then A, B, C and D are unique up to permutation and scaling of columns, meaning that any other quadruple A, B, C and D that gives rise to X is related to A, B, C, D via where is a permutation matrix, 1 , 2 , 3 and 4 are diagonal scaling matrices satisfying 1 2 3 4 = I M .

Parameter matrix estimation using QALS
In this subsection, QALS algorithm is applied to uniquely identify the four parameter matrices ) and B(f d ) for multidimensional parameter estimation. We write the observed matrix Y in (12) as follows: Based on definition 2, it is clear that (14) corresponds to one of the matrix representations of quadrilinear decomposition of a four-way array Y ∈ C M×N×2×L , with P denoting the common dimension. For clear analysis, we express the 2MN × L matrix Y as Y (1) . Thus, quadrilinear decomposition can obtain an equivalent LMN × 2 matrix representation as the following form: and an equivalent 2LN × M matrix representation where Y (1) , Y (2) , Y (3) and Y (4) denote the slice sets of the four-way array Y along the snapshot, polarization, reception and transmission ways, respectively. From definition 3 and theorem I in section 3.1, it is obvious that the k-ranks of A t , A r , G and B are k A t = k A r = P, k G = min(2, P), k B = min(L, P), respectively. If multiple targets and enough snapshots are taken into account, i.e. L ≥ P ≥ 2, then the following inequality is always satisfied: (18) which means that the four-way array Y satisfies the Kruskal-rank Theorem. Thus, A t , A r , G and B can be recovered uniquely up to permutation and scaling ambiguity.
Based on the four matrix representations Y (1) , Y (2) , Y (3) and Y (4) , the parameter matrices , η) and B(f d ) can be estimated using quadrilinear alternating least square (QALS) to fit PARAFAC model. The basic idea behind QALS is to update one parameter matrix using the least squares algorithm, conditioned on previously obtained estimates for the remaining parameter matrices that define the decomposition. This process is repeated until convergence in the least squares fit. The detailed iterative algorithm is as follows: where F denote the Frobenius norm of its matrix argument. iii) Update the k th iteration solutionĜ k (θ r , φ r , γ , η) based acquiring the least squares solution of the k th where ε is an error threshold.

Multidimensional parameter estimation
Upon obtaining the estimates ofÂ t (θ t , φ t ),Â r (θ r , φ r ), G(θ r , φ r , γ , η) andB(f d ), the following problem is to estimate multidimensional parameters of P targets from these parameter matrices. The algorithm is based on available Vandermonde structure [26] and Kronecker product structure in the matrix blocks, also uses the polarizationsensitive array processing technology [27].
Define ς p as the ratio between the second and first elements of the polarization vector g(θ rp , φ rp , γ p , η p ) in (32) ς p sin γ p cos θ rp sin φ rp e jη p + cos γ p cos φ rp sin γ p cos θ rp cos φ rp e jη p − cos γ p sin φ rp .
The estimate of ς p can be calculated bŷ whereς p is the estimate of ς p . Rewriting ς p in (33) yields ς p = tan γ p cos θ rp sin φ rp e jη p + cos φ rp tan γ p cos θ rp cos φ rp e jη p − sin φ rp .

Doppler frequency estimation
The estimate of Doppler frequency is derived from the matrixB(f d ). Define ω p e j2πf dp T s , then Since [ ω 1 , · · · , ω P ] is the ratio between the (l + 1)th and lth rows of the matrix B(f d ) for l = 1, · · · , L, perform the average of ratios yieldinĝ whereω p is the estimate of ω p . Therefore, the Doppler frequency can be estimated bŷ Upon the above analysis, multidimensional parameter θ t , φ t , θ r , φ r , γ , η and f d for the pth target can be effectively calculated, p = 1, · · · , P.

Discussion
Note that the proposed method can obtain unique parametric identification even when multiple targets have close 2D-DODs or close 2D-DOAs while they are diversely polarized. Specifically, we consider the situation that P targets have the same 2D-DOD or the same 2D-DOA, while the other parameters are different. In this case, the inequality in (18) becomes From the Kruskal-rank Theorem in section 3.1, the sufficient condition is derived herein to guarantee that the quadrilinear decomposition model is unique. Therefore, with polarization diversity, the 2D-DODs and 2D-DOAs of multiple closely spaced targets in MIMO radar can be uniquely identified, and high-resolution estimation can be achieved. It is mentioned that, based on the signal model described in (14), the trilinear decomposition-based algorithm with TALS can also be applied to identify multidimensional parameters using three matrices A t , A r G and B as the slice sets of three-way array Y. However, for trilinear decomposition, the Kruskal-rank condition can not be satisfied in the situation that P targets have the same 2D-DOD, since the inequality k A t + k A r G + k B ≥ 2P + 2 is needed to uniquely identify the three parameter matrices. http://asp.eurasipjournals.com/content/2013/1/133 Besides, robust iterative fitting of multilinear approaches such as linear programming (LP) or weighted median filtering iteration [28] can also be applied here to further yield a better solution for robust estimation in non-Gaussian noise.
The following assumption is made in the article: (i) The antennas are assumed to be ideal, identical isotropic scatterers. In practical MIMO radar systems, the effect of mutual coupling and array self-calibration methods should be considered. (ii) The target amplitude is assumed to remain constant during snapshot collection (Swerling I model). The condition can be relaxed to Swerling II model [29] when Doppler shift is known. (iii) The number of targets is assumed to be known. In practical situation, the detection of the number of targets, such as the minimum description length (MDL) criterion, should be investigated to obtain the rank of the four-way array Y. (iv) The targets are assumed to be point sources. In practical MIMO radar implementation, distributed target models and methods can be further considered.

Simulation
The simulations are conducted to verify the effectiveness of the proposed method in this section. Consider a polarimetric bistatic MIMO radar system with M 1 = M 2 = 3 and N 1 = N 2 = 3. The inter-sensor spacings are d t = d r = λ/2. The duration of a snapshot is T s = 50μs, and the length of transmit codes is K = 1024. The amplitudes of P targets are β 1 = β 2 = · · · β P = 1. The number of snapshots is L = 100, and Monte Carlo trials are 100. Simulation 1. Consider P = 3 targets with 2D-DOD, 2D-DOA, polarization parameters and Doppler frequency respectively being (θ t1 , θ t2 , 000, 3, 000, 5, 000) Hz.
First, we present the simulation results of 2D-DOD, 2D-DOA, polarizations and Doppler frequency estimation for three targets when signal-to-noise ratio (SNR) = 10 dB. The average of 100 Monte Carlo trials is shown in Table 1. From Table 1, it is observed that the estimates of seven parameters are close to their true values. The effectiveness and accuracy of the proposed algorithm can be verified from the simulation.
Further, the performance of the proposed algorithm is evaluated through root-mean-square error (RMSE). For estimating a parameter α, the RMSE of the pth target can be calculated by where N t is the number of Monte Carlo trials. Figure 2a,b,c,d,e,f,g respectively show the RMSE versus SNR for joint estimation of transmit elevation, transmit azimuth, receive elevation, receive azimuth, two polarization parameters and Doppler frequency. From Figure 2, we observe that the performance of multidimensional parameter estimation has been improved gradually with the increase of SNR. The seven parameters for each target can be efficiently obtained without peak searching, and the parameters for multiple targets can be paired automatically without wrong pairing. The additional pairing computation is eliminated using our algorithm. Based on the estimated multidimensional parameters, accurate identification and 3D localization of multiple targets can be achieved in MIMO radar.
Secondly, we assume that the two targets have the same 60°). The average result of 100 Monte Carlo trials for 2D-DOD and 2D-DOA estimation is shown in Table 3.
Thirdly, we assume that the two targets have close 2D-DODs, i.e.  Table 4.
The simulations of the above scenarios indicate that the proposed algorithm can distinguish multiple targets having close 2D-DODs and 2D-DOAs, or even the same 2D-DOD or 2D-DOA by polarization diversity. Therefore, it can uniquely identify the 2D transmit/receive angles of multiple closely spaced targets with high-resolution. However, traditional methods only contain angle parameter estimation without polarization information, leading to their performance degradation in locating multiple closely spaced targets.

Simulation 3.
In this example, the proposed QALSbased algorithm is compared with TALS-based algorithm and ESPRIT-based algorithm to evaluate the performance of 2D-DOD and 2D-DOA estimation. For the other two algorithms, we extend the conventional PARAFAC TALSbased algorithm in [19] and the ESPRIT algorithm in [5] from 1D-estimation version to 2D-estimation version in order to compare them. Similarly, a uniform rectangular transmit array and a cross-dipoled uniform rectangular receive array are used in the two algorithms to estimate multidimensional parameters based on the signal model of polarimetric MIMO radar given in (12). In the extended ALS-based algorithm, the observed matrix Y in (12) is written as , it is shown that Y is a PARAFAC trilinear decomposition model. Thus, TALS-based algorithm in [19] can be applied to recover the three parameter matrices A t (θ t , φ t ), A r (θ r , φ r , γ , η) and B(f d ). Following this, 2D-DOD, 2D-DOA as well as polarization and Doppler frequency can be estimated exploiting the Vandermonde structure and Kronecker matrix products in the three parameter matrices. The extended ESPRIT-based algorithm employs the shift invariance property of transmit and receive arrays, in which both the steering vectors of transmit and receive arrays have the form of kronecker product for the rectangular arrays, as shown in (6). The signal subspace can be obtained by eigen-value decomposition of covariance matrix of the observed matrix Y in (12). Then, multidimensional parameters can be estimated by properly partitioning the signal subspace matrix into different submatrices [16]. The procedure of pairing 2D-DODs and 2D-DOAs for multiple targets is performed similarly to [30].
In simulation 3, we consider a scenario that there are P = 2 targets. Their polarization parameters and Doppler frequency are (γ 1 , γ 2 ) = (π/10, π/4), (η 1 , η 2 ) = (π/5, 4π/5) and (f d1 , f d2 ) = (1, 000, 3, 000) Hz. We firstly evaluate these three algorithms in the condition of targets with different 2D-DODs and 2D-DOAs. Assume that they  80°). The performance comparison of transmit elevation/azimuth (2D-DOD) and receive elevation/azimuth (2D-DOA) among the three algorithms is respectively shown in Figure 3a,b,c,d, where the RMSE is the average of the two targets. Figure 3 reveals that the three algorithms have close estimation performance in this condition, with QALS-based algorithm being a little better than the other two algorithms. Next, we evaluate the three algorithms in the condition of targets with the same 2D-DOA. Assume that they are located at (θ t1 , θ t2 ) = (10°, 20°), 80°). The performance comparison is respectively shown in Figure 4a,b,c,d. Figure 4 demonstrates that all the three algorithms can identify two targets with the same 2D-DOA since the receive array manifolds for the two targets are distinguishable by polarization diversity. However, the performance of QALS-based algorithm outperforms that of the other two algorithms, especially for 2D-DOA estimation. We finally evaluate the three algorithms in the condition of targets with the same 2D-DOD. Assume that they are located at 80°). The performance comparison is respectively shown in Figure 5a,b,c,d. From Figure 5, we can see that in this condition, TALS-based algorithm and ESPRIT-based algorithm cannot work in receive elevation/azimuth estimation because of the rank deficiency. Also, for transmit elevation/azimuth estimation, the RMSE of QALS-based algorithm is smaller than that of the other two algorithms. The reason is that for the QALS-based method, the Kruskal-rank condition of qudrilinear decomposition in (42) is well satisfied in the scenario that P targets have the same 2D-DOD. Thus, multidimensional parameters can be uniquely identified, and the proposed algorithm still works well in this scenario, whereas the TALS-based algorithm and ESPRITbased algorithm have severe performance deterioration. In addition, for the ESPRIT-based algorithm, an additional multidimensional parameter pairing for multiple targets is needed.

Conclusions
In this article, we investigate multidimensional target parameter estimation in a polarimetric bistatic MIMO radar system using PARAFAC quadrilinear decomposition. The signal model has been developed and a novel algorithm for joint estimation of 2D-DOA, 2D-DOD, polarization parameters and Doppler frequency has been presented based on QALS. The algorithm has many merits: (i) it requires neither multidimensional spectral peak searching nor covariance matrix estimation and several eigen-decompositions that may bring error accumulation, which enhances the accuracy of estimation; (ii) multidimensional parameters can be well paired automatically, which reduces the complexity of additional pairing; (iii) it can distinguish multiple targets having close 2D-DODs and 2D-DOAs or even the same 2D-DOD or 2D-DOA by exploiting polarization diversity and the uniqueness of quadrilinear decomposition. Thus, the performance of multi-target 2D-DOD and 2D-DOA estimation in polarimetric MIMO radar has been greatly improved; (iv) unlike the previous algorithms with only respect to 1D-DOD and 1D-DOA estimation in bistatic MIMO radar, the proposed algorithm can obtain 2D-DOA and 2D-DOD estimation with high-resolution, which is of importance for accurate identification and 3D localization of multiple targets in MIMO radar. http://asp.eurasipjournals.com/content/2013/1/133