A Fast PARAFAC Algorithm for Parameter Estimation in Monostatic FDA-MIMO Radar

This paper studies the joint range and angle estimation of monostatic frequency diverse array multiple-input multiple-output (FDA-MIMO) radar and proposes a joint estimation algorithm. First, the transmit direction matrix is converted into real values by unitary transformation, and the Vandermonde-like matrix structure is used to construct an augmented output that doubles the aperture of the receive array. Then the augmented output is combined into a third-order tensor. Next, the factor matrices are initially estimated. Finally, the direction matrices are estimated utilizing parallel factor (PARAFAC) decomposition, and the range and angle are calculated by employing least square fitting. As contrasted with the classic PARAFAC method, the proposed method can estimate more targets and provide better estimation performance, and requires less computational complexity. The availability and excellence of the proposed method are reflected by numerical simulations and complexity analysis.


Introduction
The use of parameter estimation to locate targets is a crucial part of radar target detection [1]. However, traditional radar struggles to accomplish the more accurate parameter estimation required by modern society. As a result, researchers used multiple-input multiple-output (MIMO) technology [2][3][4][5], which has had significant success in the world of communication, in the field of radar, resulting in the formation of a MIMO radar. Due to the crossed signal transmission and reception, MIMO radar generates many virtual elements in the array, which enhances degrees of freedom and parameter estimation performance significantly [6,7]. Accordingly, target localization using parameter estimation of MIMO radar has become a hot research field in academia [8]. However, since the beam pattern of MIMO radar is only influenced as a result of the angle information, two targets with unequal ranges and equal angles cannot be identified [9].
Frequency diversity array (FDA) radar was proposed by Antonik and Wicks in 2006 [10]. By adjusting the carrier frequency of each antenna [11,12], the FDA radar controls the beam forms of various range-angle relationships, providing the possibility of joint range and angle calculation for the target. For the purpose of giving MIMO radar the range-dependent beam advantages of FDA radar, researchers proposed the FDA-MIMO radar system [13,14].
Joint range-angle parameter estimation has become a popular study direction with the employment of FDA-MIMO radar. In [15], the rotational invariance technique was used to estimate parameter values for monostatic FDA-MIMO by Li B. et al. In [16], Feilong Liu proposed an improved method for estimation of signal parameters that uses rotational invariance techniques (ESPRIT); it employs unitary transformation technology to increase the range and angle calculation accuracy of FDA-MIMO radar and decrease the complexity. However, for the joint range and angle parameter calculation, the improved ESPRIT method still suffers from the low accuracy problem. In [17], Xiong J. et al. proposed a method for parameter estimation of monostatic FDA-MIMO radar based on multiple signal classification (MUSIC). This method improves the estimation accuracy to a considerable level, but the complexity is too high to ensure efficiency in practical applications, as the algorithm involves a two-dimensional spectral peak search. The above algorithms are all algorithms based on matrix analysis, which ignores the multi-dimensional information of multi-dimensional data. To make better use of this information, researchers began to study the methods of multi-dimensional data processing.
For multi-dimensional data processing, parallel factor (PARAFAC) is a method that preserves the multidimensional structure of the data and performs well in parameter estimation [18][19][20][21][22]. Thus, scholars began to use the PARAFAC method to estimate the parameter values of FDA-MIMO radar [23][24][25]. The traditional PARAFAC algorithm randomly sets the initial values, e.g., initializing the factor matrices with Gaussian random matrices. However, setting the initial values in this way has great uncertainty and often leads to a high number of iterations for convergence, which makes the algorithm's complexity increase greatly, especially in the cases of large array numbers or high sample numbers. Consequently, how to reduce the complexity of the PARAFAC decomposition process has also become an important part of using PARAFAC. One way to decrease the complexity of PARAFAC is to employ the tensor compression technique to get a smaller tensor [26][27][28]. However, there exists information loss in the compression process, which causes degradation of performance. Therefore, reducing the complexity of the PARAFAC decomposition process by taking a suitable initial value also becomes a way to apply PARAFAC decomposition efficiently [29].
As the parameter estimation performance has a great relationship with the number of snapshots and the array aperture [30][31][32][33][34], in an attempt to increase the parameter estimation performance, many scholars began to expand the number of snapshots and the array virtual aperture of the received signal model by constructing the augmented output, so as to increase the accuracy of the parameter estimation. To better exploit the advantages of FDA-MIMO radar, the parameter calculation accuracy of the PARAFAC-like method can be improved by constructing augmented output. In this paper, we firstly propose improved fast parallel factor decomposition for joint range and angle value calculation of monostatic FDA-MIMO radar, which utilizes a structure of the signal matrix similar to that of Vandermonde to construct the augmented output, which expands the aperture of the receive array performance and improves the parameter estimation. The expansion of the receive array aperture leads to an increase in algorithm complexity. On that account, we use the idea of the propagator method (PM) [35] to obtain the initial estimated values of the factor matrices to lower the complexity of the process of the PARAFAC decomposition. Numerical simulations reflect that the proposed method is more excellent than the classic PARAFAC algorithm [23]. Utilizing the PM to reduce computational complexity does not decrease estimation performance.
For this paper, the outline is as follows. In Section 2, we introduce some tensor operations and the FDA-MIMO radar model for the parameter calculation. In Section 3, we construct the augmented output and perform the initial estimation of the factor matrices, and use parallel decomposition to estimate the angle and range. In Section 4, we report simulation experiments and performance analysis. Section 5 is the conclusion of this work.
For the notation of this paper, Table 1 gives the definitions.

Notations Definitions
identity matrix the vector outer product Khatri-Rao product D n (B) the diagonal matrix consists of the n-th row of B · F the Frobenius norm the real part operator angle(A) the phase angle of A

Tensor Operations
First, we describe some operations on third-order tensors related to this paper [36], taking X ∈ I 1 × I 2 × I 3 as an example. Operation 1: Rank-one tensor.
X is rank one if it can be formulated as: where c ∈ I 1 × 1, b ∈ I 2 × 1 and a ∈ I 3 × 1. Each element of X is the multiplication of the corresponding element of the vector: Operation 2: PARAFAC decomposition. The PARAFAC decomposition decomposes X up to the sum of the rank-one tensors.
where c k , b k and a k are factor matrices. The elements of X are:

Signal Model
In this paper, all estimation work is based on monostatic FDA-MIMO radar. Assume the MIMO radar consists of M antennas in the transmitting array, and the receiving array has N antennas. Both the transmitter and receiver are designed as uniform linear arrays (ULA), and the spacing in the transmitting array d t and the receiving array d r is equal to one-half the wavelength. Due to the limitation of sampling technology and the physical characteristics of FDA, the receiving array cannot be frequency-controlled, so generally only the transmitting array adopts a frequency-controlled array. The construction of monostatic FDA-MIMO radar is given in Figure 1. It can be found that in the monostatic FDA-MIMO radar, the direction of arrival (DOA) is identical to the direction of departure (DOD). Design the carrier frequencyf 1 of the first element as the reference frequency. Then for the m-th element, the carrier frequency is: where ∆ f represents the frequency increment between adjacent array elements and ∆ f f 1 .
The baseband signals generated by different elements of the transmitting array are orthogonal to each other; then the emission waveform of the element of the m-th element satisfies: where φ m (t) denotes the baseband waveform, T represents the duration of the radar pulse, and Q = √ E/M, where E denotes the overall transmitted energy. Assume the baseband signals are all positively orthogonal to each other. Set the time delay to T, and then normalize them these signals: Suppose there are K uncorrelated narrowband stationary targets in the far-field. The range of targets far exceeds the array aperture of the monostatic FDA-MIMO radar. The range and angle of the k-th target are given by r k and θ k , respectively. r k has to be confined to the interval (0, c/2∆ f ), where c represents the speed of electromagnetic wave propagation, due to the constraint of the maximum uncertainty range [37].
The transmitting array is an FDA array, so the transmitting steering vector contains angle information and range information. This steering vector is [38]: r(r k ) contains range information expressed as: Since the receiving array is not frequency control, the receiving steering vector only contains the angle information. This steering vector is [39]: Next, the received data of the receiver are [40]: where f s , f k and β k denote the pulse repeat frequency, Doppler frequency and the radar cross section (RCS), respectively. Since the signals send by the transmitter have orthogonality, the received data are subjected to matched filtering processing, respectively. The output result of the m-th matching filter is [40]:x where a tm (θ k , r k ) is the m-th element of a t (θ k , r k ). Stack the outputs of all matched filters into a vector , which can be expressed as [41]: Assume that a t (θ k , r k ), a r (θ k ) and the RCS are consistent for J snapshots, and S = [s(1), · · · , s(J)]. The received data are expressed as: S can be expressed as: Each row of the signal matrix S has a structure similar to the Vandermonde matrix, as can be seen.

Construction of Augmented Output
The section obtains an augmented output whose receiving array aperture is twice that of the original signal, which improves the maximum number of identifiable targets [34] and the parameter estimation performance.
Let U M be the unitary matrix. When M is odd, the unitary matrix is [42]: where M = 2k + 1.
When M is even, the unitary matrix is: where M = 2k. The transmitting steering vector can instead be rewritten as: Assuming that M is odd, U M a t (θ k , r k ) is: Define According to Equations (20) and (21), the direction steering vector is transformed by: Define the direction matrix as . According to Equation (22), the direction matrix can be transformed by: The received data can be transformed by: According to Equation (16), the following formula is obtained: where The augmented output can be constructed as [34]: where a e (θ k ) = [u k a T r (θ k ), α k u * k a H r (θ k )] T . According to Equation (26), the receive direction matrix of the augmented output is: The transmit direction matrix of the augmented output is: According to Equations (27) and (28), Equation (26) can be modified as: where X n is Under the influence of noise, the model of the augmented output is: According to the trilinear model [39], the tensor of the augmented output is: Depending on the matrix expansion of the tensor, the matrix Y is [43]: The matrix Z is [43]:

Initial Estimation of Factor Matrix
The PARAFAC technique usually employs the trilinear alternating least squares (TALS) method to directly decompose the obtained model. However, the traditional TALS algorithm directly uses random initial values to iterate, which often leads to slow convergence of the iterative decomposition process, which increases the complexity of the algorithm. The receiving array aperture was expanded by a factor of two in the previous section, considerably increasing the complexity of each iteration in the trilinear decomposition process. To decrease the complexity of trilinear decomposition, we used the PM idea for initial estimation of factor matrices, with scale variation to get better initial values, reduce the number of parallel factor iterations and decrease the complexity in the trilinear decomposition process.
After elementary row transformation, the reconstructed received augmented output X E can be expressed as: where G ∈ C 2N M×2N M is the row transformation matrix, which can be expressed as: Define V = A F A E . It is divided into two parts: where V 1K ∈ C K×K , V 2K ∈ C (2MN−K)×K , P ∈ C (2MN−K)×K and I K is the K × K identity matrix. Let P c be: Then the following relationship can be obtained: For the reconstructed received augmented output X E , the covariance matrix is: The matrix R is partitioned by: where R 1 ∈ C 2MN×K ; R 2 ∈ C 2MN×(2MN−K) . Without noise, R 1 and R 2 satisfy: In application, the covariance matrix R is: The estimate of the matrix P is the minimum of the following equation: The estimate of P is: According to Equation (38), the estimate ofP c can be obtained: Then, the matrixP c is partitioned by: According to P c V 1K = V and V = A F A E , the following formulas are obtained: There are scale ambiguity and permutation ambiguity in the PARAFAC decomposition, namely, can be expressed as: Equation (48) can be reformulated as: Equation (49) can be reformulated as: According to Equations (52) and (53), the following formula is obtained: By transforming Equation (54), we get: The eigenvalue decomposition of P † 1 P 2 is performed to obtain the estimated D 2 (Â F ) of D 2 (Ã F ) and the estimatedV 1K of V 1K . By Equation (39), the estimatedV =P cV 1K of V can be obtained. Next,V is partitioned by: where the matrices V 1 , V 2 , · · · , V M are all of the size 2N × K. From Equation (56), we can get the initial estimateÂ Eini = V 1 ofÃ E , and we can also get can be obtained: Owing toÃ F = A F D 1 (A F ), the first row ofÃ F are all 1. We can obtain the initial estimateÂ Fini ofÃ F .
According to

Trilinear Decomposition
The receive and transmit direction matrices are estimated by employing the TALS method in this Section [18].
Depending on Equation (29), the least square (LS) fitting is: The LS of matrix S is estimated as: whereÂ E represents the former estimates of A E andÂ F represents the former estimates of A F . Depending on Equation (33), the LS fitting is: The LS of matrix A E is estimated as: whereÂ F represents the former estimates of A F andŜ represents the former estimates of S.
Depending on Equation (34), the LS fitting is: The LS of matrix A E is estimated as: whereÂ E represents the former estimates of A E andŜ represents the former estimates of S.

Angle and Range Estimation
By TALS method, the final transmit direction matrixÂ F and receive direction matrix A E with scale ambiguity and permutation ambiguity can be obtained. Permutation ambiguity does not affect angle and range estimation; normalization can be used to resolve scale ambiguity.
After resolving the scale ambiguity, the range and angle estimated values can be obtained by the LS fitting. a Ek and a Fk stand for k-th column of A E and A F , respectively. We can get: where k = angle(α k (u * k ) 2 ), and normal stands for normalization. Divide A E into two parts to fit: whereĥ Ek1 andĥ Ek2 represent the first and last N rows of h Ek . Define For v Ek1 , the LS solution is as follows: For v Ek2 , the LS solution is as follows: The estimated value of the angle θ k is: The LS fit to A F is given by: Define Forv tk , the LS solution is as follows: The estimated value of the range r k is: The following are the major steps in the proposed method: Step 1: According to Equations (24) and (26), obtain the augmented matrix X E ; Step 2: According to Equation (35), obtain the reconstructed received augmented output X E ; Step 3: Get the covariance matrixR via Equation (43); Step 4: The propagation operatorP is obtained by Equation (45); Step 5: Estimate initial estimateÂ Eini with scale variation via Equations (55) and (56); Step 6: Estimate initial estimateÂ Fini with scale variation via Equations (56) and (57); Step 7: Estimate initial estimateŜ ini with scale variation via Equation (58); Step 8: According to Equation (32), obtain the third-order tensor X ; Step 9: According to Equations (60), (62) and (64), obtain the direction matrices estimateŝ A E andÂ F by iterative update; Step 10: Calculate angle via Equation (73); Step 11: Calculate range via Equation (77);

CRB
Depending on [17], the Cramér-Rao Bound(CRB) for monostatic FDA-MIMO angle and range estimation can be calculated by: where P 1 where W θ and W r are: where a k = a r (θ k ) ⊗ a t (θ k , r k ), k = 1, · · · , K.

Experimental Results and Performance Analysis
In this section, three incoherent targets are assumed to be in the far field of the monostatic FDA-MIMO radar systems. The numbers of elements of the transmitting array and receiving array are M = 7 and N = 5, respectively. The angle and range of the three targets are: (θ 1 , wherer k,l andθ k,l are the estimated values of the range and angle of k-th target for l-th Monte Carlo trial, respectively. In the first experiment, we set the number of snapshots to J = 100 and the signal-tonoise ratio (SNR) to SNR = 20. The estimated values of angle and range are close to the set values, as shown in Figure 2. The stability and availability of the proposed method are reflected in this experiment. In the second experiment, we investigate the range and angle parameter estimation performance of the proposed method under the different SNRs. The number of snapshots is set to J = 100. From Figures 3 and 4, it can be found that the proposed method outdoes the classic PARAFAC method [23], the ESPRIT method [15] and the unitary ES-PRIT method [16], and the estimation performance becomes better as SNR increases. The reason for this is that the proposed method doubles the aperture of the array. Simultaneously, the proposed method performs similarly to the improved Gaussian random initialization method.  In the third experiment, we investigate the range and angle parameter estimation performance of the proposed method under the different numbers of snapshots. The SNR is set to 20. In Figures 5 and 6, it can be found that the estimated range and angle values of the proposed method are more accurate than those of the classic PARAFAC method [23], the ESPRIT method [15] and the unitary ESPRIT method [16], and the RMSE gradually lowers as the number of snapshots rises.
In the fourth experiment, we investigate the estimation performance of the proposed algorithm for the closed-target condition case. The three closed targets are (θ 1 , r 1 ) = (30 • ,60 km), (θ 2 , r 2 ) = (35 • ,60 km) and (r 2 , θ 2 ) = (40 • ,60 km), and J is 100. In practice, if the targets are closed in space, the correlation between the column vectors of the orientation matrix becomes larger, which can make target separation difficult. In Figures 7 and 8, it can be observed that the proposed method has smaller RMSE values than the classic PARAFAC method [23], ESPRIT method [15] and unitary ESPRIT method [16] in the case of closedtarget condition.

Complexity Analysis
For the proposed method, the complexity of constructing the augmented output process is O{N 2 M 2 + MN J 2 + M 2 N 2 J}, the complexity of the initial estimation process of the factor matrix is O{8N 2 M 2 J + 4M 2 N 2 K + 6NK 2 + 2MK 3 + 8MNK 2 + 2N MK J}, the computational complexity of each iteration process in the trilinear decomposition is O{3K 2 + 6KMN J + 2K 2 (2N J + MJ + 2MN)} and the number of iterations is n 1 . The entire complexity of the proposed method is O{N 2 M 2 + MN J 2 + M 2 N 2 J + 8N 2 M 2 J + 4M 2 N 2 K + 6NK 2 + 2MK 3 + 8MNK 2 + 2N MK J + n 1 (6KMN J + 3K 3 + 2K 2 (2N J + MJ + 2MN))}. The number of iterations in the proposed method is n 1 . The overall complexity of the improved Gaussian random initialization algorithm is O{N 2 M 2 + MN J 2 + M 2 N 2 J + n 2 (6KMN J + 3K 3 + 2K 2 (2N J + MJ + 2MN))}. n 2 represents the number of iterations of the improved Gaussian random initialization algorithm. Usually, n 2 is larger than n 1 , and often n 2 is much larger than n 1 . Therefore, the overall complexity of the proposed method is less than that of the improved Gaussian random initialization method. Thus, as to verify the total complexity of the two algorithms, we performed 20,000 experiments and plotted the average time spent on each estimate by the two algorithms for these 20,000 experiments. In this experiment, the numbers of elements of the transmitter and receiver were M = 7 and N = 5, respectively. The SNR and the number of targets were 10 and 3, respectively. Figure 9 shows that the proposed method took less time than the improved Gaussian random initialization algorithm, and that the difference became more noticeable as the total number of snapshots rose. This proves that the method of estimating the factor matrix using the PM idea reduces the number of iterations of the trilinear decomposition process, thereby reducing the complexity of the algorithm.

Conclusions
In this paper, a joint range and angle estimation method based on fast improved PARAFAC decomposition for Monostatic FDA-MIMO radar was proposed. The proposed method constructs an augmented output with twice the aperture of the receive array via data expansion technique and unitary transformation. Then, the trilinear decomposition process is accelerated by factor matrix initialization and the direction matrix is obtained. Lastly, the range and angle estimated values are calculated by the LS fitting. In comparison with the classic PARARAC method, the proposed method can calculate more precise parameters and detect more targets. In addition, the proposed method uses the PM idea to calculate the factor matrix, avoiding the excessive number of iterations of the trilinear decomposition process.