Tensor-Based Reduced-Dimension MUSIC Method for Parameter Estimation in Monostatic FDA-MIMO Radar

Abstract

Frequency diverse array (FDA) radar has attracted much attention due to the angle and range dependence of the beam pattern. Multiple-input-multiple-output (MIMO) radar has high degrees of freedom (DOF) and spatial resolution. The FDA-MIMO radar, a hybrid of FDA and MIMO radar, can be used for target parameter estimation. This paper investigates a tensor-based reduced-dimension multiple signal classification (MUSIC) method, which is used for target parameter estimation in the FDA-MIMO radar. The existing subspace methods deteriorate quickly in performance with small samples and a low signal-to-noise ratio (SNR). To deal with the deterioration difficulty, the sparse estimation method is then proposed. However, the sparse algorithm has high computation complexity and poor stability, making it difficult to apply in practice. Therefore, we use tensor to capture the multi-dimensional structure of the received signal, which can optimize the effectiveness and stability of parameter estimation, reduce computation complexity and overcome performance degradation in small samples or low SNR simultaneously. In our work, we first obtain the tensor-based subspace by the high-order-singular value decomposition (HOSVD) and establish a two-dimensional spectrum function. Then the Lagrange multiplier method is applied to realize a one-dimensional spectrum function, estimate the direction of arrival (DOA) and reduce computation complexity. The transmitting steering vector is obtained by the partial derivative of the Lagrange function, and automatic pairing of target parameters is then realized. Finally, the range can be obtained by using the least square method to process the phase of transmitting steering vector. Method analysis and simulation results prove the superiority and reliability of the proposed method.

1. Introduction

MIMO radar has received widespread attention in the field of target parameter estimation, which has great development potential [1,2,3,4,5]. With constantly developing and progressing in the field of wireless communication, MIMO radar has also gained popularity in the fields of national defense, navigation, remote sensing, and unmanned driving [6,7,8,9]. It has been confirmed in practical applications that MIMO radar transmits orthogonal waveforms to form virtual array elements to enlarge the array aperture, which has a higher degree of freedom and spatial resolution [10,11,12]. Based on such advantages, MIMO radar has a high application value in overcoming the fading effect and optimizing the performance of the target location [13,14,15]. Among the wide applications of MIMO radar, target location plays an important role in both industry and research fields. For example [12], multiple unmanned aerial vehicles are used to achieve cooperative localization of maritime targets. The paper [16] proposed a DOA and polarization estimation in MIMO radar.

With the harsh electromagnetic environment, only the angle information cannot satisfy the increasing demands for target localization. Some scholars consider the range estimation as another important parameter to describe the target location more accurately. However, the beam pointing is independent of range, which makes it difficult to obtain range information with traditional MIMO radar [3]. In response to the disadvantage, FDA radar has been proposed in [17,18]. The FDA radar sets the frequency increment on the transmitting array elements so that the beam pointing changes with the target range [19]. This can be applied for the joint DOA and range estimation, which has great potential in practical applications [1,20,21].

However, for FDA-MIMO radar, the angle information and range information of the target are coupled and the phase period is ambiguous. It is possible that one target corresponds to more than one angle [22]. So as to deal with these problems, a double-pulse target location method has proposed in [23]. Firstly, the frequency increment can be preset to zero, i.e., the traditional phased array radar can be employed to estimate the angle of the target. Then the range information of targets is obtained. Moreover, the method eliminates the phase ambiguity by adjusting the frequency increment and obtains the correct target parameter. In [24,25], the subarray-based FDA has been proposed, where the transmitting array is composed of two sub-arrays, with different frequency increments of each sub-array. The angle and range coupling problem is solved by applying the nonlinear frequency offset method. In fact, FDA-MIMO radar is the most practical method for decoupling the angle and range of the target, which combines the advantages of FDA radar and MIMO radar.

Based on the idea, the [26] has proposed the FDA-MIMO radar. It has adopted the characteristics of the freedom of transmission waveform and beam range pointing, which improves the speed resolution and target observation time of the radar system [27,28]. These characteristics are very suitable for suppressing clutter interference. Currently, there are some target location methods based on FDA-MIMO radar [29,30]. A joint angle and range of target minimum variance distortion-free response method is proposed in [31], but this method cannot fully use the high DOF. In [32], the author has proposed the estimation of signal parameters via the rotational invariance technique (ESPRIT) method. To optimize the computational efficiency, the unitary ESPRIT method has been proposed in [33]. A target location method via real-valued subspace decomposition has been proposed in [34], which implements the unitary ESPRIT method in sub-array FDA-MIMO. However, these methods rely on the rotation invariant structure of the received signal and are only suitable when the frequency of the transmitting antenna increases linearly. Moreover, the subspace method relies on the accuracy of the signal subspace, while the effectiveness declines rapidly in the case of small samples or low SNR. As the singular value decomposition loses the array aperture, the performance is lower than the MUSIC method [35]. In [8], it has proposed a two-dimensional MUSIC method in FDA-MIMO radar. The method is difficult to be employed in practical applications due to the high computation complexity, which is brought by two-dimensional spectrum peak search. In addition, sparse DOA and range estimation methods attract the attention of scholars. The block sparse representation with mutual coupling has been proposed in [36]. In [37], an off-grid method based on Bayesian learning has been proposed, which implements the estimation of DOA and range. Unfortunately, high computation complexity is also a disadvantage of the sparse method, and it is hard to bear the cost of actual applications.

However, these methods introduced above have two main shortages: (1) the received signals are stacked in matrix form and the target parameters are estimated by matrix decomposition, which ignores the inherent multi-dimensional structure of the received signal [38,39]. (2) the peak search methods (MUSIC) have huge computation complexity and are not suitable to apply in practical scenarios if the work needs to deal with multi-dimensional peak search, e.g., joint DOA and range estimation. Therefore, it is necessary to study a target parameter estimation method that can effectively utilize multi-dimensional structure and overcome the performance degradation of low SNR and snapshot number. Besides, the new method also needs to reduce the dimensionality of the conventional two-dimensional MUSIC method.

In this paper, a tensor-based dimension reduction MUSIC method in FDA-MIMO radar is proposed. We use tensor storage to receive signals to capture multi-dimensional structures. Then HOSVD is used to process the signal model based on tensor and reduce the dimensionality of the MUSIC method [40]. Thus, DOA information can be estimated by a one-dimensional spectrum peak search. After that, according to the characteristic of the transmit-receive steering vector in FDA-MIMO radar, the Lagrange multiplier method is applied to estimate the transmitting steering vector. Finally, the phase of the transmitting steering vector is extracted to obtain the range of the target. As each DOA corresponds to a unique range parameter, the automatic matching work is also realized. We list the main contributions of the proposed method as follows:

(1) In the tensor domain, the proposed method realizes joint angle and range estimation of FDA-MIMO radar. The inherent multi-dimensional structure of the signal model is stored in the tensor, which improves the accuracy of target parameter estimation. This is the unique advantage of the tensor-based method.

(2) With the Lagrange multiplier method, the DOA and range parameters are decoupled, which further reduces the dimension of peak search. Then we use the one-dimensional spectrum function to estimate the DOA of the target and improve the accuracy.

(3) As the range of the target is not limited to a small interval, it’s too difficult to find the range information in a defined search area. Instead of peak search, our work uses constraint conditions to derive the transmitting steering vector. Then the phase is extracted and the least square method can be applied to estimate the range information that also reduces computation complexity. Moreover, the proposed method achieves automatic matching of the angle and range.

The definition notations are shown in

Table 1. Related notation.

Notation	Definition
$\mathcal{Q}$(bold Euler script letter)	tensor
$\mathit{A}$(bold capital letter)	matrix
$\mathit{a}$(bold lowercase letter)	vector
∘	Hadamard product
⊗	Kronecker product
⊙	Khatri-Rao product
${\mathit{I}}_Q$	$Q\times Q$ identity matrix
${\mathit{0}}_Q$	$Q\times Q$ zero matrix
${(·)}^{*}$	conjugate of matrix
${(·)}^{\mathrm{T}}$	transpose of matrix
${(·)}^{\mathrm{H}}$	conjugation-transpose of matrix
$\mathrm{diag}(\ast )$	diagonalization of matrix
$\mathrm{angle}(\ast )$	extract phase
${\mathbb{C}}^{M\times N}$	$M\times N$ matrix set

to facilitate the derivation of subsequent formulas.

2. Basic Knowledge of Tensor and Signal Model Based on Tensor

2.1. Basic Knowledge of Tensor

This section mainly defines the basic operations of tensor to facilitate the derivation and display of the proposed method. Please refer to [41,42,43] for more details.

Definition 1.

(Tensor unfolding of n-mode): Define$\mathcal{Q}\in {\mathbb{Z}}^{I_1\times I_2\cdots \times I_N}$as an N-order tensor, whose dimension is $I_1\times I_2\cdots \times I_N$. The n-mode tensor unfolding of $\mathcal{Q}$ can be expressed as ${\left[\mathcal{Q}\right]}_{\left(n\right)}$. The $(i_n,j)$-th element of ${\left[\mathcal{Q}\right]}_{\left(n\right)}$ is equal to $\left(i_1,i_2,\cdots ,i_N\right)$-th element of $\mathcal{Q}$ that $j=1+{\sum }_{k=1,k\ne n}^N(i_k-1)P_k$, and $P_k={\prod }_{m=1,m\ne n}^{k-1}{I}_m$.

Definition 2.

(Tensor-matrix product of n-mode): It can be expressed as$\mathcal{X}=\mathcal{Q}{\times }_n\mathit{R}$for the n-mode product, which is the tensor $\mathcal{Q}\in {\mathbb{Z}}^{I_1\times I_2\cdots \times I_N}$ and the matrix $\mathit{R}\in {\mathbb{Z}}^{J_n\times I_n}$. Further we derive the element correspondence between $\mathcal{X}$ and $\mathcal{Q}$ as ${\left[\mathcal{X}\right]}_{i_1,i_2,\cdots ,i_{n-1},j_n,i_{n+1},\cdots ,i_N}={\sum }_{i_n=1}^{I_n}{\left[\mathcal{Q}\right]}_{i_1,i_2,\cdots ,i_{n-1},i_n,i_{n+1},\cdots ,i_N}{\left[\mathit{R}\right]}_{j_n,i_n}$. In addition, the dimension of $\mathcal{X}$ is $I_1\times I_2\times \cdots \times I_{n-1}\times J_n\times I_{n+1}\times \cdots \times I_N$.

Definition 3.

(HOSVD of Tensor) The high-order singular value decomposition (HOSVD) is an important part of tensor decomposition [44]. It is written as

$$\mathcal{Q}=\mathcal{G}{\times }_1{\mathit{F}}_1{\times }_2{\mathit{F}}_2{\times }_3\cdots {\times }_N{\mathit{F}}_N$$

(1)

where the core tensor of $\mathcal{Q}$ is defined as $\mathcal{G}$, whose dimension is the same as to $I_1\times I_2\times \cdots \times I_N$. Then select the left singular vector of ${\left[\mathcal{Q}\right]}_{\left(n\right)}$ to form ${\mathit{F}}_n\in {\mathbb{C}}^{I_n\times I_n}(n=1,2,\cdots ,N)$.

2.2. Signal Model Based on Tensor

We present a monostatic FDA-MIMO radar system, as shown in Figure 1, which is wide beam and narrow-band. The uniform linear arrays(ULA) is applied to the transmitting array and the receiving array, and d represents the spacing between adjacent antennas of the array. Generally, we define the reference frequency as $f_1$, which is the transmitting frequency of the first element in the transmitting array. In addition, the radar system has M transmitting elements and N receiving elements. Then $f_m$ denotes the transmitting frequency of m-th transmitting element, which can be defined as

$$f_m=f_1+(m-1)\Delta f,(m=1,2,\cdots ,M)$$

(2)

where $\Delta f$ stands for the frequency increment. It should be noted that $\Delta f$ should be much smaller than $f_1$.

The transmitted signal of the m-th transmitting element is written as

$$s_m\left(t\right)={\varphi }_m\left(t\right)e^{j2\pi f_{m}t},0\le t\le T$$

(3)

where ${\varphi }_m\left(t\right)$ and T denote baseband signal of the m-th transmitting antenna and the radar pulse duration, respectively. The transmitted signal is wide beam and narrow-band complex. On the premise that the waveforms are orthogonal to each other, and the time shift is set to $\tau$ that is derived as

$${\int }_0^{\mathrm{T}}{\varphi }_m\left(t\right){\varphi }_n^{*}(t-\tau )e^{j2\pi \Delta f(m-n)t}dt=\left\{\begin{array}{c}0,m\ne n,\forall \tau \hfill \\ 1,m=n,\tau =0.\hfill \end{array}\right.$$

(4)

Suppose there are P targets, which are unrelated to each other and in the far-field. For the p-th $(p=1,2,\cdots ,P)$ target, the DOA and range parameters are denoted by ${\theta }_p$ and $r_p$, respectively. The $a_t({\theta }_p,r_p)$ represents the transmitting steering vector that is written as

$${\mathit{a}}_t({\theta }_p,r_p)=\left[\begin{array}{c}1\\ e^{-j2\pi \left(f_{1}dsin\left({\theta }_p\right)-2\Delta f{r}_p\right)/c}\\ ⋮\\ e^{-j2\pi (M-1)\left(f_{1}dsin\left({\theta }_p\right)-2\Delta f{r}_p\right)/c}\end{array}\right]$$

(5)

The $a_r\left({\theta }_p\right)$ represents the receiving steering vector that is written as

$${\mathit{a}}_r\left({\theta }_p\right)=\left[\begin{array}{c}1\\ e^{-j2\pi f_{1}dsin\left({\theta }_p\right)/c}\\ ⋮\\ e^{-j2\pi (N-1)f_{1}dsin\left({\theta }_p\right)/c}\end{array}\right].$$

(6)

The expression of received signal is as follows:

$$\mathit{x}\left(t\right)=[{\mathit{a}}_r\left({\theta }_1\right)\otimes {\mathit{a}}_t({\theta }_1,r_1),\cdots ,{\mathit{a}}_r\left({\theta }_P\right)\otimes {\mathit{a}}_t({\theta }_P,r_P)]\mathit{s}\left(t\right)+\mathit{n}\left(t\right)$$

(7)

where the speed of light is denoted by c, and $c=3\times {10}^8$ m/s.

Then arrange the received signal into the form of matrix multiplication, which is written as

$$\begin{array}{c}\hfill \mathit{x}\left(t\right)=\mathit{A}\mathit{s}\left(t\right)+\mathit{n}\left(t\right)\end{array}$$

(8)

where the transmitting signal is denoted by $\mathbf{s}\left(t\right)={[s_1\left(t\right),s_2\left(t\right),\cdots ,s_p\left(t\right)]}^{\mathrm{T}}\in {\mathbb{C}}^{P\times 1}$, whose elements are the product of the reflection coefficient of P targets and the doppler frequency shift, and the element can be written as $s_p\left(t\right)={\beta }_p\left(t\right)e^{-j2\pi f_1{r}_p/c}$. In addition, we can get $\mathit{A}={\mathit{A}}_r\odot {\mathit{A}}_t$ where ${\mathit{A}}_r=[{\mathit{a}}_r\left({\theta }_1\right),\cdots ,{\mathit{a}}_r\left({\theta }_P\right)]\in {\mathbb{C}}^{N\times P}$ and ${\mathit{A}}_t=[{\mathit{a}}_t({\theta }_1,r_1),\cdots ,{\mathit{a}}_t({\theta }_P,r_P)]\in {\mathbb{C}}^{M\times P}$ represents the direction matrix of receiving and transmitting, respectively; $\mathit{n}\left(t\right)\in {\mathbb{C}}^{NM\times 1}$ represents the received additional white Gaussian noise matrix.

The received signal matrix is formed after receiving L snapshots that is denoted by

$$\mathit{X}=[{\mathit{A}}_r\odot {\mathit{A}}_t]{\mathit{S}}^{\mathrm{T}}+\mathit{N}$$

(9)

where the specific composition of the received signal is $\mathit{X}=[\mathit{x}\left(t_1\right),\mathit{x}\left(t_2\right),\cdots ,\mathit{x}\left(t_L\right)]$, and $\mathit{S}={[\mathit{s}\left(t_1\right),\mathit{s}\left(t_2\right),\cdots ,\mathit{s}\left(t_L\right)]}^{\mathrm{T}}\in {\mathbb{C}}^{L\times P}$. After receiving L snapshots, the received noise matrix can be denoted by $\mathit{N}=[\mathit{n}\left(t_1\right),\mathit{n}\left(t_2\right),\cdots ,\mathit{n}\left(t_L\right)]\in {\mathbb{C}}^{NM\times L}$, which is uniform Gaussian white noise.

Then to store the received signal, we employ a third-order tensor $\mathcal{X}$, whose dimension is $M\times N\times L$. According to the Definition 1, Equation (9) is divided into blocks along the snapshot direction. Let ${\left[\mathcal{N}\right]}_{\left(3\right)}^{\mathrm{T}}=\mathit{N}$, then stack the blocks of $\mathit{X}$ following the third-dimension of $\mathcal{X}$, which is written as

$${\left[\mathcal{X}\right]}_{\left(3\right)}^{\mathrm{T}}=\mathit{X}.$$

(10)

In this way, we get the signal model based on tensor for FDA-MIMO radar. Using tensor algebra for target location can optimize the performance of target parameters estimation.

3. Doa and Range Estimation VIA Tensor for FDA-Mimo Radar

3.1. Signal Subspace Estimation VIA HOSVD

The tensor-based target location method requires HOSVD of the received signal, which is different from the traditional matrix decomposition-based method [43]. According to Definition 2 and Definition 3, the HOSVD of tensor-based received signal model $\mathcal{X}$ is written as

$$\mathcal{X}=\mathcal{G}{\times }_1{\mathit{U}}_1{\times }_2{\mathit{U}}_2{\times }_3{\mathit{U}}_3$$

(11)

where the core tensor can be denoted by $\mathcal{G}$, which is obtained by HOSVD. ${\mathit{U}}_n\in {\mathbb{C}}^{I_n\times I_n}$$(n=1,2,3)$ are formed of the left singular, which is obtained by the n-mode tensor unfolding of $\mathcal{X}$ according to Definition 1, namely ${\left[\mathcal{X}\right]}_{\left(n\right)}={\mathit{U}}_n{\mathrm{\Lambda }}_n{\mathit{V}}_n^{\mathrm{H}}$. Moreover, because we have assumed P targets, the rank of $\mathcal{X}$ is P. The subspace based on the tensor is obtained by truncated HOSVD, which is written as

$${\mathcal{X}}_t={\mathcal{G}}_t{\times }_1{\mathit{U}}_{t1}{\times }_2{\mathit{U}}_{t2}$$

(12)

where ${\mathcal{G}}_t=\mathcal{X}{\times }_1{\mathit{U}}_{t1}^{\mathrm{H}}{\times }_2{\mathit{U}}_{t2}^{\mathrm{H}}{\times }_3{\mathit{U}}_{t3}^{\mathrm{H}}$ denotes the truncated core tensor. Perform eigenvalue decomposition on ${\mathit{U}}_n(n=1,2,3)$ and extract the eigenvectors corresponding to the P largest values of ${\mathit{U}}_n(n=1,2,3)$ to form ${\mathit{U}}_{tn}$. Then substituting ${\mathcal{G}}_t$ into Equation (12), after simplification, the result is written as

$${\mathcal{X}}_t=\mathcal{X}{\times }_1\left({\mathit{U}}_{t1}{\mathit{U}}_{t1}^{\mathrm{H}}\right){\times }_2\left({\mathit{U}}_{t2}{\mathit{U}}_{t2}^{\mathrm{H}}\right){\times }_3{\mathit{U}}_{t3}^{\mathrm{H}}.$$

(13)

Therefore, the subspace based on tensor that is expressed as

$${\mathit{U}}_s={\left[{\mathcal{X}}_t\right]}_{\left(3\right)}^{\mathrm{T}}=({\mathit{U}}_{t2}{\mathit{U}}_{t2}^{\mathrm{H}}\otimes {\mathit{U}}_{t1}{\mathit{U}}_{t1}^{\mathrm{H}}){\left[\mathcal{X}\right]}_{\left(3\right)}^{\mathrm{T}}{\mathit{U}}_{t3}^{*}$$

(14)

where ${\left[\mathcal{X}\right]}_{\left(3\right)}={\mathit{U}}_3{\mathrm{\Lambda }}_3{\mathit{V}}_3^{\mathrm{H}}$. Similarly, ${\left[\mathcal{X}\right]}_{\left(3\right)}^{\mathrm{T}}\approx {\mathit{V}}_{t3}^{*}{\mathrm{\Lambda }}_{t3}{\mathit{U}}_{t3}^{\mathrm{T}}$ is obvious. Substituting ${\left[\mathcal{X}\right]}_{\left(3\right)}^{\mathrm{T}}$ into Equation (14), simplifying the formula is written as [38]

$${\mathit{U}}_s=({\mathit{U}}_{t2}{\mathit{U}}_{t2}^{\mathrm{H}}\otimes {\mathit{U}}_{t1}{\mathit{U}}_{t1}^{\mathrm{H}}){\mathit{V}}_{t3}^{*}{\mathrm{\Lambda }}_{t3}.$$

(15)

So far, the signal subspace ${\mathit{U}}_s$ based on tensor has been estimated.

3.2. DOA Estimation VIA Tensor-Based Reduced-Dimension Music

The signal subspace and the noise subspace are orthogonal to each other, which is the basis of the MUSIC method [45]. The signal subspace is transformed into a unit orthogonal column matrix by orthogonalization. After orthogonal projection, we can get the noise subspace, which can be expressed as

$${\mathit{U}}_{noise}{\mathit{U}}_{noise}^{\mathrm{H}}={\mathit{I}}_{NM}-{\mathit{U}}_o{{\mathit{U}}_o}^{\mathrm{H}}$$

(16)

The tensor-based spectrum function is constructed as

$$f(\theta ,r)=\frac{1}{{[{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{a}}_t(\theta ,r)]}^{\mathrm{H}}[{\mathit{I}}_{NM}-{\mathit{U}}_o{{\mathit{U}}_o}^{\mathrm{H}}][{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{a}}_t(\theta ,r)]}$$

(17)

where ${\mathit{U}}_o$ represents the orthogonal basis of ${\mathit{U}}_s$, and define ${\mathit{U}}_{orth}={\mathit{U}}_{noise}{\mathit{U}}_{noise}^{\mathrm{H}}$. However, Equation (17) involves peak searching for two-dimensional spectral peaks, which requires a lot of computation and is inefficient in practical application. Inspired by the Ref. [40], a dimensional strategy is investigated to reduce the computation complexity based on the characteristic of transmit-receive steering vector of FDA-MIMO radar, and the detail is showed as follows.

In order to facilitate the derivation of the subsequent formula, according to the properties of the Kronecker product [43]. Set ${\mathit{W}}_1\in {\mathbb{C}}^{d_1\times d_2}$, ${\mathit{W}}_2\in {\mathbb{C}}^{d_2\times d_3}$, ${\mathit{W}}_3\in {\mathbb{C}}^{d_4\times d_5}$ and ${\mathit{W}}_4\in {\mathbb{C}}^{d_5\times d_6}$, we can get

$$\begin{array}{cc}\hfill \left({\mathit{W}}_1{\mathit{W}}_2\right)\otimes \left({\mathit{W}}_3{\mathit{W}}_4\right)& =({\mathit{W}}_1\otimes {\mathit{W}}_3)({\mathit{W}}_2\otimes {\mathit{W}}_4)\hfill \end{array}$$

(18)

$$\begin{array}{c}\hfill {\mathit{W}}_1\otimes {\mathit{W}}_3=\left({\mathit{W}}_1{\mathit{I}}_{d_2}\right)\otimes \left({\mathit{I}}_{d_4}{\mathit{W}}_3\right)=({\mathit{W}}_1\otimes {\mathit{I}}_{d_4})({\mathit{I}}_{d_2}\otimes {\mathit{W}}_3)\end{array}$$

(19)

Because of ${\mathit{a}}_r\left(\theta \right)\in {\mathbb{C}}^{N\times 1}$ and ${\mathit{a}}_t(\theta ,r)\in {\mathbb{C}}^{M\times 1}$, we can get the simplified ${\mathit{a}}_r\left(\theta \right)\otimes {\mathit{a}}_t(\theta ,r)$, which is expressed as

$$\begin{array}{cc}\hfill {\mathit{a}}_r\left(\theta \right)\otimes {\mathit{a}}_t(\theta ,r)& =\left[{\mathit{a}}_r\left(\theta \right){\mathit{I}}_1\right]\otimes \left[{\mathit{I}}_M{\mathit{a}}_t(\theta ,r)\right]\hfill \\ \hfill & =[{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{I}}_M]\left[{\mathit{I}}_1{\mathit{a}}_t(\theta ,r)\right]\hfill \\ \hfill & =[{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{I}}_M]{\mathit{a}}_t(\theta ,r).\hfill \end{array}$$

(20)

So we define and simplify as follows

$$\begin{array}{cc}\hfill \mathit{F}(\theta ,r)& ={[{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{a}}_t(\theta ,r)]}^{\mathrm{H}}{\mathit{U}}_{orth}[{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{a}}_t(\theta ,r)]\hfill \\ \hfill & ={\mathit{a}}_t{(\theta ,r)}^{\mathrm{H}}{[{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{I}}_M]}^{\mathrm{H}}{\mathit{U}}_{orth}[{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{I}}_M]{\mathit{a}}_t(\theta ,r)\hfill \\ \hfill & ={\mathit{a}}_t{(\theta ,r)}^{\mathrm{H}}\mathit{F}\left(\theta \right){\mathit{a}}_t(\theta ,r)\hfill \end{array}$$

(21)

For the convenience of expression, we define $\mathit{F}\left(\theta \right)={[{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{I}}_M]}^{\mathrm{H}}{\mathit{U}}_{orth}[{\mathit{a}}_r\left(\theta \right)\otimes {\mathit{I}}_M]$. Inspired by the Ref. [40], the quadratic optimization of Equation (21) when constraint condition is ${\mathit{e}}_o^{\mathrm{H}}{\mathit{a}}_t(\theta ,r)=1$ and ${\mathit{e}}_o={[1,0,\cdots ,0]}^{\mathrm{T}}\in {\mathbb{C}}^{M\times 1}$, so as to eliminate the solution of ${\mathit{a}}_t(\theta ,r)={\mathit{0}}_M$. Because Equation (21) is essentially an issue of conditional extreme value. We introduce the Lagrange multiplier method to address this issue, which can be constructed as follows

$$L(\theta ,r)={\mathit{a}}_t{(\theta ,r)}^{\mathrm{H}}\mathit{F}\left(\theta \right){\mathit{a}}_t(\theta ,r)+\lambda ({\mathit{e}}_o^{\mathrm{H}}{\mathit{a}}_t(\theta ,r)-1)$$

(22)

where $\lambda$ is a constant.

According to the nature of Lagrangian extreme value, the partial derivative of the Equation (22), which can be expressed as

$$\begin{array}{cc}\hfill \frac{\partial L(\theta ,r)}{\partial {\mathit{a}}_t(\theta ,r)}& =2\mathit{F}\left(\theta \right){\mathit{a}}_t(\theta ,r)-\lambda {\mathit{e}}_o=0\hfill \\ \hfill \Rightarrow {\mathit{a}}_t(\theta ,r)& =\frac{\lambda }{2}\mathit{F}{\left(\theta \right)}^{-1}{\mathit{e}}_o\hfill \end{array}$$

(23)

According to the constraint condition, we can get

$$\begin{array}{cc}\hfill & {\mathit{e}}_o^{\mathrm{H}}{\mathit{a}}_t(\theta ,r)=1\hfill \\ \hfill \Rightarrow & {\mathit{a}}_t(\theta ,r)=\frac{1}{{\mathit{e}}_o^{\mathrm{H}}}\hfill \end{array}$$

(24)

Substituting Equation (24) into Equation (23), it is shown as

$$\begin{array}{cc}\hfill & \frac{\lambda }{2}\mathit{F}{\left(\theta \right)}^{-1}{\mathit{e}}_o={\mathit{a}}_t(\theta ,r)=\frac{1}{{\mathit{e}}_o^{\mathrm{H}}}\hfill \\ \hfill \Rightarrow & \lambda =\frac{2}{{\mathit{e}}_o^{\mathrm{H}}\mathit{F}{\left(\theta \right)}^{-1}{\mathit{e}}_o}\hfill \end{array}$$

(25)

So ${\mathit{a}}_t(\theta ,r)$ is derived as

$${\mathit{a}}_t(\theta ,r)=\frac{\mathit{F}{\left(\theta \right)}^{-1}{\mathit{e}}_o}{{\mathit{e}}_o^{\mathrm{H}}\mathit{F}{\left(\theta \right)}^{-1}{\mathit{e}}_o}.$$

(26)

Inserting ${\mathit{a}}_t(\theta ,r)$ of Equation (26) into Equation (17) and Equation (21), we can estimate the DOA via

$$\begin{array}{cc}\hfill \widehat{\theta }& =\arg \max f(\theta ,r)\hfill \\ & =\arg \min {\mathit{a}}_t{(\theta ,r)}^{\mathrm{H}}\mathit{F}\left(\theta \right){\mathit{a}}_t(\theta ,r)\hfill \\ & =\arg \min {\mathit{e}}_o^{-1}\mathit{F}\left(\theta \right){\mathit{e}}_o^{-\mathrm{H}}\hfill \\ & =\arg \max {\mathit{e}}_o^{\mathrm{H}}\mathit{F}{\left(\theta \right)}^{-1}{\mathit{e}}_o\hfill \end{array}$$

(27)

It is clearly concluded from Equation (27) that the spectrum function has actually been independent of the range, thus achieving dimension reduction. DOAs of the P targets should correspond to the P largest peaks of $f(\theta ,r)$, which is obtained by $\theta \in [-{90}^{\circ },{90}^{\circ }]$.

3.3. Range Estimation

Since the DOA and range are coupled in ${\mathit{a}}_t(\theta ,r)$, if the conventional MUSIC method is used, the computation complexity is too high due to the excessive range. We use reduced-dimension to estimate the range of the target.

Defining the estimation of DOA as ${\widehat{\theta }}_p(p=1,2,\cdots ,P)$, and then substituting DOA into Equation (26) to estimate the corresponding transmitting steering vector ${\widehat{\mathit{a}}}_t({\theta }_p,r_p)$. Normalize ${\widehat{\mathit{a}}}_t({\theta }_p,r_p)$ to eliminate scale ambiguity. Extract the phase of ${\widehat{\mathit{a}}}_t({\theta }_p,r_p)$ and combine it with Equation (5), which can be written as

$$\begin{array}{cc}\hfill {\widehat{\mathit{\varphi }}}_p& =-\mathrm{angle}\left(\widehat{\mathit{a}}t({\theta }_p,r_p)\right)\hfill \\ \hfill & =\left[\begin{array}{c}0\\ \pi (sin\left({\theta }_p\right)-4\Delta f\frac{r_p}{c})\\ ⋮\\ \pi (M-1)(sin\left({\theta }_p\right)-4\Delta f\frac{r_p}{c})\end{array}\right].\hfill \end{array}$$

(28)

obtain the range information of the target, we introduce the least square method. The least-square is constructed as follows

$$\min \left|\right|\mathit{G}{\mathit{q}}_p-{\widehat{\mathit{\varphi }}}_p{\left|\right|}_F^2$$

(29)

where ${\mathit{q}}_p\in {\mathbb{C}}^{2\times 1}$ represents the estimation parameter vector. $\mathit{G}$ can be expressed as

$$\mathit{G}={\left[\begin{array}{cccc}1& 1& \cdots & 1\\ 0& \pi & \cdots & (M-1)\pi \end{array}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{M\times 2}$$

(30)

Define ${\widehat{\mathit{q}}}_p$ as the estimation of ${\mathit{q}}_p$, which can be expressed as

$${\widehat{\mathit{q}}}_p={\left({\mathit{G}}^{\mathrm{T}}\mathit{G}\right)}^{-1}{\mathit{G}}^{\mathrm{T}}{\widehat{\mathit{\varphi }}}_p$$

(31)

Define the second element of ${\mathit{q}}_p$ as ${\widehat{\mathit{q}}}_p\left(2\right)$, which is the estimated value of $sin\left({\theta }_p\right)-4\Delta f\frac{r_p}{c}$. Substituting the DOA ${\widehat{\theta }}_p$ into Equation (31) to obtain the range information of target that can be written as

$${\widehat{r}}_p=\frac{c(sin\left({\widehat{\theta }}_p\right)-{\widehat{\mathit{q}}}_p\left(2\right))}{4\Delta f}.$$

(32)

So far, using the proposed method, we have estimated the DOA and range parameters of the target.

4. Performance Analysis of the Proposed Method

4.1. Computation Complexity

We introduce the computation complexity to reflect the efficiency of the proposed method, which can be expressed as follows:

• (1) The HOSVD computation complexity of $\mathcal{X}\in {\mathbb{C}}^{M\times N\times L}$ is $O\left(MNL\right(M+N+L\left)\right)$ in Equation (11);
• (2) The signal subspace estimation needs $O\left(4PLMN\right)$ in Equation (15);
• (3) In Equation (21), the dimensionality reduction of the two-dimensional search requires $O\left(M^2{N}^2{P}^2(MN+P^2)\right)$;
• (4) The spectrum peak search of DOA estimation in Equation (27) is $O(d_{\theta }\left(MP\right)!(MP-1))$, where $d_{\theta }$ represents the search times within the search DOA, and $(*)!$ stands for factorial;
• (5) Computing the range requires $O(2M^{3}P+M^{2}P)$;

The total computation complexity is $O\{4MNL+M^2{N}^2{P}^2(MN+P^2)+d\left(MP\right)!(MP-1)+2M^{3}P+M^{2}P+MNL(M+N+L)\}$.

We choose ESPRIT [32], MUSIC, [8] and Tensor-ESPRIT [30] as comparison methods. According to Refs. [32,46] and combined with calculations, we can get the computation complexity of the ESPRIT, MUSIC, and Tensor-ESPRIT.

Table 2. Computation complexity comparison.

Method	Computation Complexity
Proposed	$O\{4MNL+M^2{N}^2{P}^2(MN+P^2)+d_{\theta }\left(MP\right)!(MP-1)$ $+2M^{3}P+M^{2}P+MNL(M+N+L)\}$
ESPRIT	$O\{{\left(2MN\right)}^{2}L+{\left(2MN\right)}^3+4(5MN-2M-2N)K^2$ $+K^3(L+M)+MN{K}^2+31K^3\}$
MUSIC	$O\{d_{\theta }*d_{range}\left[{\left(MN\right)}^2(2(MN-P)+MN)\right]$ $+{\left(MN\right)}^{2}L+M^{2}N+M{N}^2\}$
Tensor-ESPRIT	$O\{2{\left(MNL\right)}^3+MNL(M+N+L)+MLK(N+K)$ $+K^3(L+M)+MN{K}^2+31K^3\}$

shows the computation complexity of the proposed method, ESPRIT, MUSIC, and Tensor-ESPRIT.

It should be noted that in Table 2, $d_{range}$ represents the search times of range-dimension. The conventional target range is kilometer level, so $d_{range}$ is very large. It causes the computation complexity of MUSIC to be not at the same level as that of other methods. Subsequent simulations confirmed this view.

4.2. Cram$\acute{\mathit{e}}$r-Rao Bound

The received signal can be converted into

$$\begin{array}{cc}\hfill \mathit{X}& =[{\mathit{A}}_r\odot {\mathit{A}}_t]{\mathit{S}}^{\mathrm{T}}+\mathit{N}=\mathit{A}{\mathit{S}}^{\mathrm{T}}+\mathit{N}.\hfill \end{array}$$

(33)

Define $\mathit{a}({\theta }_p,r_p)={\mathit{a}}_r\left({\theta }_p\right)\otimes {\mathit{a}}_t({\theta }_p,r_p)$. Let $\zeta$ represent the signal power, and the noise power is denoted by ${\sigma }^2$.

The corresponding Fisher information matrix (FIM) is expressed as

$${\mathit{F}}_{\theta }=2L\zeta Re\left[\left({\mathit{Q}}_{\theta }^{\mathrm{H}}{\mathit{W}}_A^{\perp }{\mathit{Q}}_{\theta }\right)\circ {\mathit{R}}_n^{-1}\right]$$

(34)

$${\mathit{F}}_r=2L\zeta Re\left[\left({\mathit{Q}}_r^{\mathrm{H}}{\mathit{W}}_A^{\perp }{\mathit{Q}}_r\right)\circ {\mathit{R}}_n^{-1}\right]$$

(35)

where ${\mathit{R}}_n={\sigma }^2{\mathit{I}}_P$ represents the noise covariance matrix, ${\mathit{Q}}_{\theta }=\left[\frac{\partial \mathit{a}({\theta }_1,r_1)}{\partial {\theta }_1},\cdots ,\frac{\partial \mathit{a}({\theta }_p,r_p)}{\partial {\theta }_p}\right]$, ${\mathit{Q}}_r=\left[\frac{\partial \mathit{a}({\theta }_1,r_1)}{\partial r_1},\cdots ,\frac{\partial \mathit{a}({\theta }_p,r_p)}{\partial r_p}\right]$ and ${\mathit{W}}_A^{\perp }={\mathit{I}}_{\mathit{MN}}-\mathit{A}{\left({\mathit{A}}^{\mathrm{H}}\mathit{A}\right)}^{-1}{\mathit{A}}^{\mathrm{H}}$.

The partial derivative can be expressed as

$$\frac{\partial \mathit{a}({\theta }_p,r_p)}{\partial {\theta }_p}=\frac{\partial {\mathit{a}}_r\left({\theta }_p\right)}{\partial {\theta }_p}\otimes {\mathit{a}}_t({\theta }_p,r_p)+a_r\left({\theta }_p\right)\otimes \frac{\partial {\mathit{a}}_t({\theta }_p,r_p)}{\partial {\theta }_p}$$

(36)

$$\frac{\partial \mathit{a}({\theta }_p,r_p)}{\partial r_p}=a_r\left({\theta }_p\right)\otimes \frac{\partial {\mathit{a}}_t({\theta }_p,r_p)}{\partial r_p}$$

(37)

with

$$\frac{\partial {\mathit{a}}_r\left({\theta }_p\right)}{\partial {\theta }_p}=-j2\pi \frac{f_1}{c}dcos\left({\theta }_p\right){\mathit{D}}_N{\mathit{a}}_r\left({\theta }_p\right)$$

(38)

$$\frac{\partial {\mathit{a}}_t({\theta }_p,r_p)}{\partial {\theta }_p}=-j2\pi \frac{f_1}{c}dcos\left({\theta }_p\right){\mathit{D}}_M{\mathit{a}}_t({\theta }_p,r_p)$$

(39)

$$\frac{\partial {\mathit{a}}_t({\theta }_p,r_p)}{\partial r_p}=-j4\pi \frac{\Delta f}{c}{\mathit{D}}_M{\mathit{a}}_t({\theta }_p,r_p)$$

(40)

where ${\mathit{D}}_K=\mathrm{diag}(0,\cdots ,K-1)$.

The ${\mathrm{CRB}}_{\theta }$ can be expressed as

$${\mathrm{CRB}}_{\theta }=\mathrm{mean}\left(\mathrm{diag}\left({\mathit{F}}_{\theta }^{-1}\right)\right)$$

(41)

$${\mathrm{CRB}}_r=\mathrm{mean}\left(\mathrm{diag}\left({\mathit{F}}_r^{-1}\right)\right)$$

(42)

where $\mathrm{mean}(\ast )$ represents the mean of the array.

5. Numerical Simulations

In this section, we perform numerical simulations to analyze the effectiveness and reliability of the proposed method under different situations. The FDA-MIMO radar related parameter setting are as follows: $M=N=8$; $f_1=10\mathrm{GHz}$; $\Delta f=1000\mathrm{Hz}$$d=\frac{c}{2f_1}$. The $P=3$ independent targets are supposed at $({\theta }_1,r_1)=(10.{923}^{\circ },58$ km), $({\theta }_2,r_2)=(50.{456}^{\circ },30$ km), $({\theta }_3,r_3)=(-5.{789}^{\circ },9$ km), respectively. The radar and target parameters remain unchanged unless otherwise stated. The total of Monte Carlo trials can be denoted by T, and $T=500$.

5.1. Spectrum Peak Search for DOA Estimation

The parameters of this simulation are as follows: the number of snapshots to $L=50$, and $SNR$ = 20 dB. Select the matrix decomposition-based classical MUSIC method as the comparison method. Comparing the performance of DOA estimation of the two methods. The spatial spectrum is shown in Figure 2.

In the actual simulation of the proposed method, we firstly set the search range and step size to $[-{90}^{\circ },{90}^{\circ }]$ and $1^{\circ }$, respectively, to get the angle information ${\mathrm{DOA}}_{1^{\circ }}$. Then marrow down the two search parameters appropriately, namely $[{\mathrm{DOA}}_{1^{\circ }}-5^{\circ },{\mathrm{DOA}}_{1^{\circ }}+5^{\circ }]$ and $0.1^{\circ }$. In this way, the computation complexity is reduced.

It can be found from Figure 2 that the accuracy of the signal subspace is improved, making the spatial spectrum more sharp and accurate. The reason is that the proposed method utilizes the multi-dimensional structure of the tensor-based received data to improve the estimation accuracy of the signal subspace. In order to show the difference of the spatial spectrum more clearly, the more high-precision spatial spectrum is shown in Figure 3. Figure 3 shows the part of the spatial spectrum at three target angles with the accuracy of $0.1^{\circ }$. It can be seen from Figure 3 that the proposed method has more significant peaks and narrower bandwidth, which also proves that the proposed method has a higher resolution.

5.2. 2D Point Cloud of the Target Landing Point

In this simulation, the estimation of results are shown in Figure 4. The signal parameters are set to $L=50$ and $SNR$ = 20 dB. The DOA and range of target are displayed as X-axis and Y-axis, respectively. We can get from the Figure 4 that the landing point of the estimated target coincides with the real target position. The automatic matching of DOA and range can be realized. This can intuitively reflect the superiority of the stability and accuracy of the proposed method.

5.3. RMSE Performance

It is important for the root mean square error (RMSE) to measure the comprehensive performance of the proposed method. The important index can be expressed as

$$RMS{E}_{\theta }=\frac{{\sum }_{p=1}^P\sqrt{{\sum }_{t=1}^T{({\widehat{\theta }}_p^t-{\theta }_p)}^2/T}}{P}$$

(43)

$$RMS{E}_r=\frac{{\sum }_{p=1}^P\sqrt{{\sum }_{t=1}^T{({\widehat{r}}_p^t-r_p)}^2/T}}{P}$$

(44)

where the estimation of DOA and range are denoted as ${\widehat{\theta }}_p^t$ and ${\widehat{r}}_p^t$, which are the t-th monte carlo tarils estimation results of the p-th target, respectively. The ESPRIT, MUSIC, Tensor-ESPRIT, and CRB are used to measure the comprehensive performance of the proposed method.

We show the performance of the methods with SNR when $L=50$. The RMSE of DOA and range are shown in Figure 5 and Figure 6, respectively. We can obtain that the proposed method outperforms all other methods in accuracy and stability. This is because the proposed method uses a multi-dimensional structure through tensor, which improves the accuracy of the target location. It’s worth noting that the ESPRIT-based method loses the array aperture when decomposing and estimating the signal subspace, so the performance is lower than that of the MUSIC-based method.

Figure 7 and Figure 8 show the RMSE of DOA and range estimation of these methods with L when $SNR$ = 5 dB. It corresponds to the DOA RMSE and range RMSE, respectively. Similarly, we introduce the comparison methods. We can get that the RMSE curve of the proposed method is smoother and closer to CRB. This shows that the performance of the proposed method is the best. The performance of the ESPRIT-MUSIC method will surpass that of MUSIC when the total of snapshots is high. It is because the increase in the number of snapshots compensates for the loss of the array aperture during signal decomposition. However, the proposed method uses tensor to ensure performance. In addition, the proposed method can also achieve excellent target positioning at low snapshots.

5.4. Probability of Successful Detection

The probability of successful detection (PSD) is another important index, which can measure the stability and robustness of the method.

Define the total of successful detections as V. PSD is denoted by

$$PSD=V/T\times 100\%$$

(45)

If the difference between the estimated and the true of DOA and range are within $\pm 0.1^{\circ }$ and ± 0.1 km, respectively, then it can be regarded as a victory target detection. It can be regarded as a successful detection, when the difference between the estimated value and the true value of DOA and range is within $\pm 0.1^{\circ }$ and $\pm 0.1$ km, respectively.

Figure 9 and Figure 10 show the PSD of DOA and range estimation with SNR when the number of snapshots is $L=50$, respectively. We can conclude that the method is the best compared to the comparison methods. This indicates that the stability and effectiveness of the proposed method is excellent. In addition, the proposed method can achieve a 100% successful detection rate with a lower SNR, which indicates that our method overcomes the defect of performance degradation when SNR is low.

5.5. The Simulation Time Versus Trial Number

In this simulation, we compare the simulation running time of the methods, and the results are shown in Figure 11. The signal parameters are set to $SNR$ = 20 dB and $L=50$. In particular, the total of Monte Carlo trials is $T=100$. The X-axis and Y-axis represent the trial number and simulation running time, respectively.

Figure 11 shows the simulation running time comparison of ESPRIT, MUSIC, Tensor-ESPRIT, and Proposed method. From the previous calculation of computation complexity, it can be seen that the computation complexity of MUSIC is much higher than that of the other methods, which is caused by the range-dimensional spectral peaks. The proposed method eliminates the high computational cost of range-dimension spectral peaks through reduced-dimension. It can be seen from the partial enlargement that the simulation running time of the proposed time is close to that of ESPRIT and Tensor-ESPRIT methods. This shows that the computation complexity of the proposed method is the same order of magnitude as that of ESPRIT and Tensor-ESPRIT methods. The performance of the proposed method is more excellent under the same computation complexity.

6. Conclusions

We propose a reduced-dimensionality MUSIC method based on tensor in this paper, which realizes joint DOA and range estimation in FDA-MIMO radar. A tensor-based signal model can be established to store inherent multi-dimensional structures that can optimize the performance and stability of target parameter estimation. Then use the HOSVD to obtain the signal subspace based on tensor, which can be used to construct the two-dimensional spectrum function. The Lagrange multiplier method is utilized to reduce the dimensionality of the spectrum function. The estimation of DOA is realized by a one-dimensional spectrum peak search. The constraint condition is constructed to eliminate the limit situation when the transmitting steering vector is zero matrix. We can derive the transmitting steering vector from the partial derivative of the Lagrange function. Finally, extract the phase of the transmitting steering vector. Range information is estimated by the least square method. Moreover, the high computation complexity of two-dimensional MUSIC is overcomes in the proposed method, and uses tensor-based signal model to optimize the performance of the proposed method, which can remain stable in low SNR and low snapshot. The superiority of the proposed method can be verified by numerical analysis and simulation. In the future, we will optimize the proposed method to make target location more accurate and suitable for more application.