Target Parameter Estimation Algorithm Based on Real-Valued HOSVD for Bistatic FDA-MIMO Radar

Abstract

Since there is a frequency offset between each adjacent antenna of FDA radar, there exists angle-range two-dimensional dependence in the transmitter. For bistatic FDA-multiple input multiple output (MIMO) radar, range-direction of departure (DOD)-direction of arrival (DOA) information is coupled in transmitting the steering vector. How to decouple the three information has become the focus of research. Aiming at the issue of target parameter estimation of bistatic FDA-MIMO radar, a real-valued parameter estimation algorithm based on high-order-singular value decomposition (HOSVD) is developed. Firstly, for decoupling DOD and range in transmitter, it is necessary to divide the transmitter into subarrays. Then, the forward–backward averaging and unitary transformation techniques are utilized to convert complex-valued data into real-valued data. The signal subspace is obtained by HOSVD, and the two-dimensional spatial spectral function is constructed. Secondly, the dimension of spatial spectrum is reduced by the Lagrange algorithm, so that it is only related to DOA, and the DOA estimation is obtained. Then the frequency increment between subarrays is used to decouple the DOD and range information, and eliminate the phase ambiguity at the same time. Finally, the DOD and range estimation automatically matched with DOA estimation are obtained. The proposed algorithm uses the multidimensional structure of high-dimensional data to promote performance. Meanwhile, the proposed real-valued tensor-based method can effectively cut down the computing time. Simulation results verify the high efficiency of the developed method.

1. Introduction

Multiple input multiple output (MIMO) radar was developed in 2004, which can make up for the drawbacks of phased array (PA) [1,2,3]. Different from PA radar, the transmitted waveforms of MIMO radar are orthogonal to each other [4,5,6]. When the receiver completes the matched filtering, it can produce a large number of virtual array elements, which can dramatically promote radar performance [7,8,9]. However, the advantages of MIMO radar in range estimation are not prominent. Therefore, in 2006, Antonik et al. proposed frequency diversity array (FDA) [10,11]. With the development of array signal processing, some scholars developed FDA-MIMO radar [12,13]. There are two categories of FDA-MIMO: statistical radar [14,15] and collocated radar [16,17]. In this paper, the bistatic FDA-MIMO in collocated radar is taken as the research object.

FDA-MIMO radar adds frequency increment in transmitter antennas, the transmitter waveform is affected by both range and angle [18,19]. The transmitting waveform has two-dimensional dependence on range and angle [20,21]. Therefore, FDA-MIMO can estimate angle and range concurrently [22,23]. Moreover, since FDA-MIMO radar adds the information of range dimension, the degree of freedom (DOF) is increased [24,25]. Thus, FDA-MIMO radar can achieve more tasks in different environments [26]. Therefore, FDA-MIMO radar can be applied to achieve parameter estimation.

Since FDA-MIMO can provide more DOFs of the system, it can provide more target parameter information. Therefore, many parameter estimation studies for FDA-MIMO have been performed. In [27], an algorithm based on the multiple signal classification (MUSIC) method to achieve target angle-range is proposed. Since the estimation method requires two spectral peak searches to obtain angle and range, the estimation method has a lot of operational redundancy. In order to cut down the operation time, a two-stage algorithm via rotation invariance technique (ESPRIT) approach is developed to obtain the angle-range estimation [28]. The angle and range estimated by this algorithm are automatically matched. In [29], in order to cut down more operation time. Liu et al. constructed the unitary matrix to transform the subspace into real-valued data, further reducing the operational redundancy. However, the precision of this approach will be declined in low signal-to-noise ratio (SNR). Therefore, a tensor-based FDA-MIMO radar algorithm is proposed [30], this algorithm realizes spatial spectrum dimensionality reduction through Lagrange multiplier method, thus reducing the computational complexity of spectral peak search part, and achieving decoupling of angle and range. Xu et al. proposed an algorithm based on high-order-singular value decomposition (HOSVD), which retains the multi-dimensional structure of data and can obtain superior result. Although the above algorithms can realize the estimation of FDA-MIMO, the accuracy of estimation will sharply decline in low SNR and snapshot. However, in bistatic FDA-MIMO radar, these algorithms will fail because of the coupling problem between the DOD and the range of the target. At present, the research of target parameter estimation for bistatic FDA-MIMO radar is still insufficient, and the existing algorithms are mostly carried out for monostatic FDA-MIMO radar.

In this paper, a real-valued parameter estimation approach based on HOSVD is developed. This algorithm solves the problem of three-dimensional direction of arrival (DOA)-direction of departure (DOD)-range estimation in bistatic FDA-MIMO radar. Firstly, for eliminating the coupling between DOD and range, the transmitter is divided into several subarrays. Then a three-dimensional tensor data model is constructed. The original tensor is converted into a real-valued tensor by unitary transformation technique. The HOSVD algorithm is employed to obtain signal subspace. The spatial spectrum function is constructed through obtained subspace. Then the one-dimensional spatial spectrum only related to DOA information is obtained by Lagrange algorithm. The DOA is estimated by one-dimensional spatial spectral. Utilizing the constructed transmitting subarray, the DOD and range information are decoupled and the periodic ambiguity of phase is eliminated. Finally, the automatically matched DOD and range are estimated. The developed approach preserves the multidimensional structure. Furthermore, the matched DOA, DOD and range estimation can be obtained through the reduced-dimension MUSIC algorithm. The presented method not only has high precision, but also cut down the computational complexity.

In summary, the contributions of the developed algorithm are summarized as:

(1)

The developed approach can achieve joint DOA-DOD-range estimation for bistatic FDA-MIMO. The tensor signal subspace is obtained by HOSVD method. The original structure of the received data is preserved, which can greatly improve the estimation accuracy;

(2)

The proposed approach is a real-valued operation, and it utilizes the reduced-dimension MUSIC algorithm, which estimates DOA by utilizing one-dimensional spatial spectrum. It greatly reduces operational redundancy while ensuring the performance advantages;

(3)

The presented method eliminates the coupling of DOD information and range information by subarray division of transmitter. Accurate DOD and range estimations are achieved.

The notations are presented in

Table 1. Notations.

Notation	Definition
${(·)}^{\ast }$	conjugate
${(·)}^T$	transpose
${(·)}^H$	conjugate-transpose
${(·)}^{-1}$	inverse
${(·)}^{†}$	pseudo-inverse
⊗	Kronecker product
⊙	Khatri–Rao product
○	outer product
${\bigsqcup }_n$	the concatenation along the n-th mode
${\mathit{I}}_Q$	$Q\times Q$ identity matrix
${\mathbf{0}}_Q$	$Q\times Q$ zero matrix
$\lfloor ·\rfloor$	floor operator
$(·)!$	factorial
$diag(·)$	diagonalization of matrix
$angle(·)$	the phase of array elements

.

2. Signal Model

The bistatic FDA-MIMO radar is taken as the research object. From Figure 1, we can know that there are M transmitting antennas and N receiving antennas, all of which are uniform linear arrays (ULA) with half wavelength spacing. $d_t$ and $d_r$ are antenna spacing. On the basis of the definition of FDA-MIMO radar, there is a frequency increment between antennas at transmitter. The carrier frequency of the m-th antenna is given by [31]

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

(1)

where $f_0$ represents carrier frequency. $\Delta f$ stands for frequency increment. The signal transmitted by m-th transmitting antenna can be defined as

$$\begin{array}{c}\hfill s_m\left(t\right)=\sqrt{\frac{E}{M}}{\psi }_m\left(t\right)e^{j2\pi f_{m}t},\\ \hfill 0\le t\le T,m=1,2,\cdots ,M,\end{array}$$

(2)

where E stands for energy, ${\psi }_m\left(t\right)$ represents transmitting waveform. T stands for delay. Since transmitting waveforms are orthogonal to each other, the following expression is given

$${\int }_T{\psi }_m\left(t\right){\psi }_n^{\ast }(t-\tau )e^{j2\pi (m-n)\Delta ft}\mathrm{d}t=\left\{\begin{array}{cc}0,\hfill & m\ne n,\forall \tau \hfill \\ 1,\hfill & m=n,\forall \tau .\hfill \end{array}\right.$$

(3)

Assume K targets in the far-field. DOA, DOD, and range of k-th target are written as ${\theta }_k$, ${\phi }_k$, and $r_k$. The data of the k-th target received by n-th receiving antenna can be defined as [19,32]

$$y_n\left(t\right)={\rho }_k\sum _{m=1}^M{\psi }_m\left(t-{\tau }_{m,k}^t-{\tau }_{n,k}^r\right)e^{j2\pi f_m(t-{\tau }_{m,k}^t-{\tau }_{n,k}^r)},$$

(4)

where ${\rho }_k$ is complex-valued reflection coefficient of the k-th target. ${\tau }_{m,k}^t$ and ${\tau }_{n,k}^r$ stand for the time delay, which can be given by

$$\begin{array}{cc}\hfill {\tau }_{m,k}^t& =\left(r_{1,k}^t-(m-1)d_{t}sin\left({\theta }_k\right)\right)/c,\hfill \\ \hfill {\tau }_{n,k}^r& =\left(r_{1,k}^r-(n-1)d_{r}sin\left({\phi }_k\right)\right)/c,\hfill \end{array}$$

(5)

where $c=3\times {10}^8$ m/s.

The received snapshot is given by [33]

$$\mathit{X}=\sum _{k=1}^K{\xi }_k{\mathit{a}}_R\left({\theta }_k\right){\mathit{a}}_T^T\left(r_k,{\phi }_k\right)+\mathit{F},$$

(6)

$${\xi }_k={\rho }_k{e}^{j2\pi f_0{r}_k/c},$$

(7)

where $\mathit{F}$ represents noise vector.

Then ${\mathit{a}}_T(r_k,{\phi }_k)$ is defined as [34]

$${\mathit{a}}_T(r_k,{\phi }_k)=\mathit{r}\left(r_k\right)\odot \mathit{d}\left({\phi }_k\right)\in {\mathbb{C}}^{M\times 1},$$

(8)

$$\begin{array}{c}\hfill \mathit{r}\left(r_k\right)={\left[1,e^{-j2\pi \Delta f{r}_k/c},\dots ,e^{-j2\pi \Delta f(M-1)r_k/c}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{M\times 1},\end{array}$$

(9)

$$\begin{array}{c}\hfill \mathit{d}\left({\phi }_k\right)={\left[1,e^{j2\pi d_t{f}_{0}sin{\phi }_k/c},\dots ,e^{j2\pi d_t{f}_0(M-1)sin{\phi }_k/c}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{M\times 1}.\end{array}$$

(10)

From Equation (8), we can know that the range and DOD information are coupled with each other. However, the DOA information of target is only related to ${\mathit{a}}_R\left({\theta }_k\right)$, which is given by

$$\begin{array}{c}\hfill {\mathit{a}}_R\left({\theta }_k\right)={\left[1,e^{j2\pi d_r{f}_{0}sin{\theta }_k/c},\dots ,e^{j2\pi d_r{f}_0(N-1)sin{\theta }_k/c}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{N\times 1}.\end{array}$$

(11)

Therefore, the DOA estimation can be obtained from ${\mathit{a}}_R\left({\theta }_k\right)$. However, to obtain DOD and range estimation, ${\mathit{a}}_T(r_k,{\phi }_k)$ needs to be decoupled. Therefore, the transmitter is converted to P subarrays. Since each subarray is independent, the frequency increment of subarrays is unequal. The subarray ${\mathit{a}}_{TS}\left(r_k,{\phi }_k\right)$ can be written as

$${\mathit{a}}_{TS}\left(r_k,{\phi }_k\right)=\left[\begin{array}{c}{\mathit{a}}_{TS}^1\left(r_k,{\phi }_k\right)\\ {\mathit{a}}_{TS}^2\left(r_k,{\phi }_k\right)\\ ⋮\\ {\mathit{a}}_{TS}^P\left(r_k,{\phi }_k\right)\end{array}\right].$$

(12)

The frequency of mth element in pth ($p=1,2,\cdots ,P$) subarray is defined as $f_{TS}^{m,p}$, which can be given by

$$f_{TS}^{m,p}=f_{TS}^{1,p}+(m-1)\Delta f_p,m=1,2,\cdots ,M_{TS}^p,$$

(13)

where $\Delta f_p$ stands for frequency increment of pth subarray.

Figure 2 is subarray model. It can be seen that the transmitter is converted to P subarrays. Therefore, the steering vector of the pth subarray is given by

$$\begin{array}{c}\hfill {\mathit{a}}_{TS}^p(r_k,{\phi }_k)=e^{j\frac{2\pi }{c}\left(\left({\displaystyle \sum _{q=1}^{p-1}}{M}_{TS}^q\right)d_t{f}_{1}sin\left({\phi }_k\right)+(f_1-f_{TS}^{1,p}{r}_k)\right)}\times \left[\begin{array}{c}1\\ e^{j\frac{2\pi }{c}\left(d_t{f}_{1}sin\left({\phi }_k\right)-\Delta f_p{r}_k\right)}\\ ⋮\\ e^{j\left(M_{TS}^p-1\right)\frac{2\pi }{c}\left(d_t{f}_{1}sin\left({\phi }_k\right)-\Delta f_p{r}_k\right)}\end{array}\right].\end{array}$$

(14)

In this paper, the number of antennas in each subarray is equal, that is, $M_{TS}^p=M/P,(p=1,2,\cdots ,P)$. In addition, the frequency of the last antenna in the former subarray is identical to that of first element in latter subarray. In other words, $f_{TS}^{M_{TS}^p,p}=f_{TS}^{1,p+1}$, where $f_{TS}^{1,1}=f_1$, $\Delta f_1<\Delta f_2\cdots <\Delta f_P$.

Equation (6) can be converted to tensor form, which can be given by

$$\begin{array}{cc}\hfill & \mathcal{X}(n,m,l)=\sum _{k=1}^K{\mathit{A}}_R(n,k)\circ {\mathit{A}}_{TS}(m,k)\circ \mathit{S}(l,k)+{\mathit{F}}_{n,m,l},\hfill \end{array}$$

(15)

where ${\mathit{A}}_T(m,k)$ is $(m,k)$-th element of ${\mathit{A}}_T$, ${\mathit{A}}_R(n,k)$ denotes $(n,k)$-th element of ${\mathit{A}}_R$, $\mathit{S}(l,k)$ represents $(l,k)$-th element of $\mathit{S}$. $\mathit{S}={[{\xi }_1,{\xi }_2,\dots ,{\xi }_k]}^T$. ${\mathit{F}}_{n,m,l}$ stands for noise tensor, while l denotes number of snapshots.

3. The Proposed Method

3.1. Real-Valued Signal Subspace Estimation

For changing $\mathcal{X}$ into a centro-Hermitian tensor, forward–backward averaging technique is used

$$\mathcal{Z}=\left[{\mathcal{X}}_{{\bigsqcup }_3}\left({\mathcal{X}}^{\ast }\times {}_1{\mathbf{\Gamma }}_N\times {}_2{\mathbf{\Gamma }}_M\times {}_3{\mathbf{\Gamma }}_L\right)\right].$$

(16)

where the element on the anti-diagonal of ${\mathbf{\Gamma }}_n$ is 1, and the other elements are 0. Through the unitary transformation, a new real-valued tensor is obtained

$${\mathcal{Z}}_{real}=\mathcal{Z}\times {}_1{\mathbf{E}}_N^{\mathrm{H}}\times {}_2{\mathbf{E}}_M^{\mathrm{H}}\times {}_3{\mathbf{E}}_{2L}^{\mathrm{H}}$$

(17)

where the unitary matrices are defined as

$${\mathbf{E}}_K=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{\mathbf{I}}_K\hfill & j{\mathbf{I}}_K\\ {\mathbf{\Gamma }}_K\hfill & -j{\mathbf{\Gamma }}_K\end{array}\right]$$

(18)

$${\mathbf{E}}_{2K+1}=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{\mathbf{I}}_K& 0& j{\mathbf{I}}_K\\ {\mathbf{0}}_{1\times K}& \sqrt{2}& {\mathbf{0}}_{1\times K}\\ {\mathbf{\Gamma }}_K& 0& -j{\mathbf{\Gamma }}_K\end{array}\right]$$

(19)

Firstly, HOSVD algorithm is employed for ${\mathcal{Z}}_{real}$, which can be written as

$${\mathcal{Z}}_{real}={\mathcal{G}}_{real}{\times }_1{\mathit{U}}_1{\times }_2{\mathit{U}}_2{\times }_3{\mathit{U}}_3,$$

(20)

where ${\mathcal{G}}_{real}\in {\mathbb{C}}^{M\times N\times L}$ is core tensor, ${\mathit{U}}_1\in {\mathbb{C}}^{M\times M}$, ${\mathit{U}}_2\in {\mathbb{C}}^{N\times N}$ and ${\mathit{U}}_3\in {\mathbb{C}}^{L\times L}$ are composed of left singular of the mode-n tensor unfolding of ${\mathcal{Z}}_{real}$, respectively. In other words, ${\left[{\mathcal{Z}}_{real}\right]}_{\left(n\right)}={\mathit{U}}_n{\mathbf{\Lambda }}_n{\mathit{V}}_n^{\mathrm{H}}(n=1,2,3)$. Then the truncated HOSVD is employed to obtain signal subspace estimation, which is given by

$${\mathcal{Z}}_s={\mathcal{G}}_s{\times }_1{\mathit{U}}_{s1}{\times }_2{\mathit{U}}_{s2},$$

(21)

where ${\mathcal{G}}_s={\mathcal{Z}}_{real}{\times }_1{\mathit{U}}_{s1}^{\mathrm{H}}{\times }_2{\mathit{U}}_{s2}^{\mathrm{H}}{\times }_3{\mathit{U}}_{s3}^{\mathrm{H}}$ is the truncated core tensor. ${\mathit{U}}_{sn}(n=1,2,3)$ is composed of the column vector of ${\mathit{U}}_n$.

Substituting ${\mathcal{G}}_s$ into Equation (21), Equation (21) is given by

$${\mathcal{Z}}_s={\mathcal{Z}}_{real}{\times }_1\left({\mathit{U}}_{s1}{\mathit{U}}_{s1}^{\mathrm{H}}\right){\times }_2\left({\mathit{U}}_{s2}{\mathit{U}}_{s2}^{\mathrm{H}}\right){\times }_3{\mathit{U}}_{s3}^{\mathrm{H}}.$$

(22)

According to [35,36], the properties of the mode product is given by

$$\left\{\begin{array}{c}\mathcal{J}\times {}_i\mathbf{W}\times {}_j\mathbf{G}=\mathcal{J}{\times }_j\mathbf{W}{\times }_i\mathbf{G},j\ne i\hfill \\ \mathcal{J}\times {}_i\mathbf{W}\times {}_i\mathbf{G}=\mathcal{J}{\times }_i\left(\mathbf{W}\mathbf{G}\right),j=i\hfill \end{array}\right.$$

(23)

$$\begin{array}{cc}\hfill & {\left[\mathcal{J}\times {}_1{\mathbf{W}}_1{\times }_2{\mathbf{W}}_2\times \cdots {\times }_I{\mathbf{W}}_I\right]}_{\left(i\right)}=\hfill \\ \hfill & {\mathbf{W}}_i·{\left[\mathcal{J}\right]}_{\left(i\right)}·{\left[{\mathbf{W}}_I\otimes \dots \otimes {\mathbf{W}}_{i+1}\otimes {\mathbf{W}}_{i-1}\otimes \dots \otimes {\mathbf{W}}_2\otimes {\mathbf{W}}_1\right]}^{\mathrm{T}}\hfill \end{array}$$

(24)

where $\mathcal{J}$ is a tensor, W and G are matrices.

According to Equations (23) and (24), the signal subspace ${\mathit{U}}_s$ is obtained by mode-3 unfolding of ${\mathcal{Z}}_s$.

$${\mathit{U}}_s={\left[{\mathcal{Z}}_s\right]}_{\left(3\right)}^{\mathrm{T}}=({\mathit{U}}_{s2}{\mathit{U}}_{s2}^{\mathrm{H}}\otimes {\mathit{U}}_{s1}{\mathit{U}}_{s1}^{\mathrm{H}}){\left[{\mathcal{Z}}_{real}\right]}_{\left(3\right)}^{\mathrm{T}}{\mathit{U}}_{s3}^{\ast },$$

(25)

where ${\left[{\mathcal{Z}}_{real}\right]}_{\left(3\right)}={\mathit{U}}_3{\mathbf{\Lambda }}_3{\mathit{V}}_3^{\mathrm{H}}$. In addition, ${\left[{\mathcal{Z}}_{real}\right]}_{\left(3\right)}^{\mathrm{T}}\approx {\mathit{V}}_{s3}^{\ast }{\mathbf{\Lambda }}_{s3}{\mathit{U}}_{s3}^{\mathrm{T}}$. Therefore, Equation (25) can be rewritten as

$${\mathit{U}}_s=({\mathit{U}}_{s2}{\mathit{U}}_{s2}^{\mathrm{H}}\otimes {\mathit{U}}_{s1}{\mathit{U}}_{s1}^{\mathrm{H}}){\mathit{V}}_{s3}^{\ast }{\mathbf{\Lambda }}_{s3}.$$

(26)

Therefore, the tensor-based signal subspace estimation has been obtained.

3.2. DOA Estimation

According to MUSIC algorithm [37], the noise subspace is achieved by orthogonal projection, which can be expressed as [38]

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

(27)

where ${\mathit{U}}_o$ is the orthogonal basis of ${\mathit{U}}_s$.

The spectrum function is given by

$$f(\theta ,\phi ,r)=\frac{1}{{[{\mathbf{E}}_N^H{\mathit{a}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H{\mathit{a}}_{TS}\left(r_k,{\phi }_k\right)]}^{\mathrm{H}}[{\mathit{I}}_{NM}-{\mathit{U}}_o{{\mathit{U}}_o}^{\mathrm{H}}][{\mathbf{E}}_N^H{\mathit{a}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H{\mathit{a}}_{TS}\left(r_k,{\phi }_k\right)]}.$$

(28)

From Equation (28), we can know that parameter estimation can be obtained by three-dimensional spectral peak searching. In order to cut down the computing redundancy, the Lagrange multiplier approach is employed to cut down the dimension of the spectrum function [39].

On the basis of the Kronecker product [40], the expression of ${\mathbf{E}}_N^H{\mathit{a}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H{\mathit{a}}_{TS}\left(r,\phi \right)$ can be simplified, which is given by

$$\begin{array}{cc}\hfill {\mathbf{E}}_N^H{\mathit{a}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H{\mathit{a}}_{TS}\left(r,\phi \right)& =\left[{\mathbf{E}}_N^H{\mathit{a}}_R\left(\theta \right){\mathit{I}}_1\right]\otimes \left[{\mathbf{E}}_M^H{\mathit{a}}_{TS}\left(r,\phi \right)\right]\hfill \\ \hfill & =\left[{\widehat{\mathit{a}}}_R\left(\theta \right){\mathit{I}}_1\right]\otimes \left[{\mathbf{E}}_M^H{\mathit{a}}_{TS}\left(r,\phi \right)\right]\hfill \\ \hfill & =[{\widehat{\mathit{a}}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H]\left[{\mathit{I}}_1{\mathit{a}}_{TS}\left(r,\phi \right)\right]\hfill \\ \hfill & =[{\widehat{\mathit{a}}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H]{\mathit{a}}_{TS}\left(r,\phi \right).\hfill \end{array}$$

(29)

where ${\widehat{\mathit{a}}}_R\left(\theta \right)={\mathbf{E}}_N^H{\mathit{a}}_R\left(\theta \right)$. Then according to Equation (28), $\mathit{F}(\theta ,\phi ,r)$ is defined as

$$\begin{array}{cc}\hfill \mathit{F}(\theta ,\phi ,r)& ={[{\mathbf{E}}_N^H{\mathit{a}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H{\mathit{a}}_{TS}\left(r,\phi \right)]}^{\mathrm{H}}{\mathit{U}}_{orth}[{\mathbf{E}}_N^H{\mathit{a}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H{\mathit{a}}_{TS}\left(r,\phi \right)]\hfill \\ \hfill & ={\mathit{a}}_{TS}{\left(r,\phi \right)}^{\mathrm{H}}{[{\widehat{\mathit{a}}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H]}^{\mathrm{H}}{\mathit{U}}_{orth}[{\widehat{\mathit{a}}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H]{\mathit{a}}_{TS}\left(r,\phi \right)\hfill \\ \hfill & ={\mathit{a}}_{TS}{\left(r,\phi \right)}^{\mathrm{H}}\mathit{F}\left(\theta \right){\mathit{a}}_{TS}\left(r,\phi \right),\hfill \end{array}$$

(30)

where $\mathit{F}\left(\theta \right)={[{\widehat{\mathit{a}}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H]}^{\mathrm{H}}{\mathit{U}}_{orth}[{\widehat{\mathit{a}}}_R\left(\theta \right)\otimes {\mathbf{E}}_M^H]$, ${\mathit{U}}_{orth}={\mathit{U}}_{noise}{\mathit{U}}_{noise}^{\mathrm{H}}$.

The quadratic optimization problem needs to be considered in Equation (30). That is, the following constraints are considered [39], it is given by

$$\begin{array}{cc}\hfill & {\mathit{e}}_o^{\mathrm{H}}{\mathit{a}}_{TS}\left(r,\phi \right)=1\hfill \\ \hfill \Rightarrow & {\mathit{a}}_{TS}\left(r,\phi \right)={\left({\mathit{e}}_o^{\mathrm{H}}\right)}^{-1},\hfill \end{array}$$

(31)

where ${\mathit{e}}_o^{\mathrm{H}}={[1,0,\cdots ,0]}^{\mathrm{T}}\in {\mathbb{C}}^{M\times 1}$.

For solving the extreme value issue in Equation (30), Lagrange multiplier approach is employed. The Lagrange function of Equation (30) is given by

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

(32)

where $\lambda$ is Lagrangian multiplier.

The partial derivative of Equation (32) is given by

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

(33)

Substituting Equation (31) into Equation (33), we can obtain

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

(34)

Therefore, ${\mathit{a}}_{TS}\left(r,\phi \right)$ is given by

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

(35)

Finally, the DOA estimation can be obtained through Equations (28) and (35).

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

(36)

The result of DOA estimation is to search the reduced dimension MUSIC spectral peak. The first K largest peaks obtained are the DOAs of corresponding K targets.

3.3. DOD and Range Estimation

By using the DOA estimation and combining with Equation (35), the estimation of the transmit steering vector ${\widehat{\mathit{a}}}_T(r_k,{\phi }_k)$ is achieved. The scale ambiguity of ${\widehat{\mathit{a}}}_T(r_k,{\phi }_k)$ should be eliminated by normalization process. Therefore, the phase of ${\widehat{\mathit{a}}}_T(r_k,{\phi }_k)$ is expressed as

$${\mathbf{\Phi }}_T^p=diag\left[\begin{array}{c}{e}^{j\frac{2\pi }{c}[d_t{f}_{1}sin{\phi }_1-\Delta f_1{r}_1]}\\ e^{j\frac{2\pi }{c}[d_t{f}_{1}sin{\phi }_2-\Delta f_2{r}_2]}\\ ⋮\\ e^{j\frac{2\pi }{c}[d_t{f}_{1}sin{\phi }_k-\Delta f_K{r}_k]}\end{array}\right].$$

(37)

The phase of the kth diagonal element in ${\mathbf{\Phi }}_T^p$ corresponds to kth target, which is

$$\left\{\begin{array}{cc}\hfill \frac{2\pi }{c}{d}_t{f}_{1}sin\left({\phi }_k\right)-\frac{2\pi }{c}\Delta f_1{r}_k& =angle\left({\mathit{\varphi }}_{TS}^{k,1}\right)-2z_1\pi \hfill \\ \hfill \frac{2\pi }{c}{d}_t{f}_{1}sin\left({\phi }_k\right)-\frac{2\pi }{c}\Delta f_2{r}_k& =angle\left({\mathit{\varphi }}_{TS}^{k,2}\right)-2z_2\pi \hfill \\ \hfill & ⋮\hfill \\ \hfill \frac{2\pi }{c}{d}_t{f}_{1}sin\left({\phi }_k\right)-\frac{2\pi }{c}\Delta f_K{r}_k& =angle\left({\mathit{\varphi }}_{TS}^{k,P}\right)-2z_P\pi ,\hfill \end{array}\right.$$

(38)

where ${\mathit{\varphi }}_{TS}^{k,p}$ corresponds to the phase of the kth target in pth subarray, $z_i\in \mathbb{Z},i=1,2,\cdots ,P$. The value of $z_i$ cannot be determined due to the phase period ambiguity of $angle\left({\mathit{\varphi }}_{TS}^{k,p}\right)$. Since $\Delta f_1<\Delta f_2<\cdots <\Delta f_P$, subtract the subformula in Equation (38), which is given by

$$\left\{\begin{array}{cc}\hfill \frac{2\pi }{c}(\Delta f_2-\Delta f_1)r_k=& angle\left({\mathit{\varphi }}_{TS}^{k,1}\right)-angle\left({\mathit{\varphi }}_{TS}^{k,2}\right)+2\pi (z_2-z_1)\hfill \\ \hfill \frac{2\pi }{c}(\Delta f_3-\Delta f_2)r_k=& angle\left({\mathit{\varphi }}_{TS}^{k,2}\right)-angle\left({\mathit{\varphi }}_{TS}^{k,3}\right)+2\pi (z_3-z_2)\hfill \\ \hfill & ⋮\hfill \\ \hfill \frac{2\pi }{c}(\Delta f_P-\Delta f_{P-1})r_k=& angle\left({\mathit{\varphi }}_{TS}^{k,P-1}\right)-angle\left({\mathit{\varphi }}_{TS}^{k,P}\right)+2\pi (z_P-z_{P-1}).\hfill \end{array}\right.$$

(39)

It can be seen from Equation (39) that only range information is included in Equation (39). However, there is still phase period ambiguity in Equation (39). In order to eliminate phase period ambiguity, we define

$${\mathit{h}}_k=\left[\begin{array}{c}{h}_{k,1}\\ h_{k,2}\\ ⋮\\ h_{k,P-1}\end{array}\right]=\left[\begin{array}{c}\frac{\Delta f_2-\Delta f_1}{c}\\ \frac{\Delta f_3-\Delta f_2}{c}\\ ⋮\\ \frac{\Delta f_P-\Delta f_{P-1}}{c}\end{array}\right],$$

(40)

$$\begin{array}{c}\hfill {\zeta }_{k,i}=angle\left({\mathit{\varphi }}_{TS}^{k,i}\right)-angle\left({\mathit{\varphi }}_{TS}^{k,i+1}\right)+(z_{i+1}-z_i),(i=1,2,...,P-1).\end{array}$$

(41)

Then Equation (39) can be rewritten as

$$\left\{\begin{array}{cc}\hfill h_{k,1}\ast & r_k={\zeta }_{k,1}\hfill \\ \hfill h_{k,2}\ast & r_k={\zeta }_{k,2}\hfill \\ \hfill & ⋮\hfill \\ \hfill h_{k,P-1}\ast & r_k={\zeta }_{k,P-1}\hfill \end{array}\right.\Rightarrow {\mathit{h}}_k{r}_k={\mathit{\zeta }}_k,$$

(42)

where ${\mathit{\zeta }}_k={[{\zeta }_{k,1},{\zeta }_{k,2},\cdots ,{\zeta }_{k,P-1}]}^{\mathrm{T}}$.

Equation (42) can be simplified as

$$\begin{array}{cc}\hfill & 2\pi \frac{\Delta f_{TS}^{i+1}-\Delta f_{TS}^i}{c}{r}_k\hfill \\ \hfill =& angle\left({\mathit{\varphi }}_{TS}^{k,i}\right)-angle\left({\mathit{\varphi }}_{TS}^{k,i+1}\right)+2\pi (z_{i+1}-z_i).\hfill \end{array}$$

(43)

For determining the range parameter, the following conditions need to be met

$$0<2\pi \frac{\Delta f_{TS}^{i+1}-\Delta f_{TS}^i}{c}{r}_k⩽2\pi .$$

(44)

Since the range $r_k$ is positive and $\Delta f_{TS}^1<\Delta f_{TS}^2\cdots <\Delta f_{TS}^P$, $2\pi \frac{\Delta f_{TS}^{i+1}-\Delta f_{TS}^i}{c}{r}_k>0$. When $2\pi \frac{\Delta f_{TS}^{i+1}-\Delta f_{TS}^i}{c}{r}_k⩽2\pi$, $(z_{i+1}-z_i)$ has a unique solution, which can be expressed as

$$z_{i+1}-z_i=\left\{\begin{array}{c}\hfill 1,angle\left({\mathit{\varphi }}_{TS}^{k,i}\right)<angle\left({\mathit{\varphi }}_{TS}^{k,i+1}\right)\\ \hfill 0,angle\left({\mathit{\varphi }}_{TS}^{k,i}\right)>angle\left({\mathit{\varphi }}_{TS}^{k,i+1}\right).\end{array}\right.$$

(45)

According to Equation (44), the range $r_k$ can be written as

$$r_k⩽\frac{c}{\Delta f_{TS}^{i+1}-\Delta f_{TS}^i}.$$

(46)

Combining Equations (42) and (46), the range $r_k$ can be expressed as

$$r_k⩽\frac{c}{max\left(\Delta f_{TS}^{i+1}-\Delta f_{TS}^i\right)}.$$

(47)

Therefore, the effective range of radar estimation is affected by radar frequency.

By substituting Equation (45) into Equation (42) and employing the least square (LS) approach, the range estimation ${\widehat{r}}_k$ can be obtained, which is given by

$${\widehat{r}}_k={\mathit{h}}_k^{†}{\mathit{\zeta }}_k,$$

(48)

where the ${\widehat{r}}_k$ is the range estimation of kth target.

By substituting ${\widehat{r}}_k$ into Equation (39), we can obtain

$$\left\{\begin{array}{cc}\hfill \frac{2\pi }{c}{d}_t{f}_{1}sin\left({\phi }_k\right)=& angle\left({\mathit{\varphi }}_{TS}^{k,1}\right)-2z_1\pi +\frac{2\pi }{c}\Delta f_{TS}^1{r}_k\hfill \\ \hfill \frac{2\pi }{c}{d}_t{f}_{1}sin\left({\phi }_k\right)=& angle\left({\mathit{\varphi }}_{TS}^{k,2}\right)-2z_2\pi +\frac{2\pi }{c}\Delta f_{TS}^2{r}_k\hfill \\ \hfill & ⋮\hfill \\ \hfill \frac{2\pi }{c}{d}_t{f}_{1}sin\left({\phi }_k\right)=& angle\left({\mathit{\varphi }}_{TS}^{k,P}\right)-2z_P\pi +\frac{2\pi }{c}\Delta f_{TS}^P{r}_k.\hfill \end{array}\right.$$

(49)

Since $2d_t{f}_1/c⩽1$, the parameter $z_i(i=1,2,\cdots ,P)$ can be obtained, which is expressed as

$$z_i=⌊\frac{angle\left({\mathit{\varphi }}_{TS}^{k,i}\right)-\frac{2\pi }{c}\Delta f_{TS}^i{\widehat{r}}_k+\pi }{2\pi }⌋,i=1,2,\cdots ,P.$$

(50)

To simplify Equation (49), we define

$${\mathit{h}}_{\phi }=\left[\begin{array}{c}\frac{2\pi }{c}{f}_1{d}_t\\ \frac{2\pi }{c}{f}_1{d}_t\\ ⋮\\ \frac{2\pi }{c}{f}_1{d}_t\end{array}\right]\in {\mathbb{C}}^{P\times 1},$$

(51)

$${\mathbf{\psi }}_{\phi }=\left[\begin{array}{c}angle\left({\mathit{\varphi }}_{TS}^{k,1}\right)-2z_1\pi +\frac{2\pi }{c}\Delta f_{TS}^1{r}_k\\ angle\left({\mathit{\varphi }}_{TS}^{k,2}\right)-2z_2\pi +\frac{2\pi }{c}\Delta f_{TS}^2{r}_k\\ ⋮\\ angle\left({\mathit{\varphi }}_{TS}^{k,P}\right)-2z_P\pi +\frac{2\pi }{c}\Delta f_{TS}^P{r}_k\end{array}\right]\in {\mathbb{C}}^{P\times 1}.$$

(52)

Then Equation (49) can be rewritten as

$${\mathit{h}}_{\phi }sin{\phi }_k={\mathbf{\psi }}_{\phi }.$$

(53)

Utilizing the total LS approach, the DOD ${\widehat{\phi }}_k$ of kth target can be estimated, which is given by

$${\widehat{\phi }}_k=arcsin\left[{\left({\mathit{h}}_{\phi }^{\mathrm{T}}{\mathit{h}}_{\phi }\right)}^{-1}{\mathit{h}}_{\phi }^{\mathrm{T}}{\mathbf{\psi }}_{\phi }\right].$$

(54)

The DOA, DOD, and range obtained by the proposed approach are automatically matched without additional matching process.

4. Algorithm Analysis

4.1. Algorithm Summary

The developed method for bistatic FDA-MIMO radar can be simplified in Algorithm 1.

Algorithm 1 Target parameter estimation algorithm based on real-valued HOSVD for bistatic FDA-MIMO radar.

1: ${\mathit{a}}_T(r_k,{\phi }_k)$ is converted to subarray by Equation (12), 2: Signal can be converted to tensor form $\mathcal{X}$, 3: Construct the real-valued tensor ${\mathcal{Z}}_{real}$ by Equation (17), 4: Calculate ${\mathit{U}}_s$ by Equation (26), 5: Construct reduced dimensional spatial spectrum function by Equation (30), 6: $\widehat{\theta }$ is achieved by Equation (36), 7: Decouple DOD and range information by Equation (39), 8: Eliminate phase period ambiguity by Equation (40), 9: ${\widehat{r}}_k$ is obtained by Equation (48), 10: ${\widehat{\phi }}_k$ is estimated by Equation (54).

4.2. Computational Complexity

In order to highlight the advantages of the developed approach, the computational complexity of our method is given by

(1)

The computational complexity of HOSVD for $\mathcal{X}\in {\mathbb{C}}^{M\times N\times L}$ is $O(1/4MNL(M+N+L\left)\right)$ in Equation (20);

(2)

The computational complexity of signal subspace estimation is $O\left(KLMN\right)$ in Equation (26);

(3)

The computational complexity of dimensionality reduction for three-dimensional spatial spectrum in Equation (30) is $O\left(M^2{N}^2{K}^2(MN+K^2)\right)$;

(4)

The computational complexity of spatial spectrum search for DOA estimation in Equation (36) is $O(d_{\theta }\left(MK\right)!(MK-1))$, where $d_{\theta }$ is the DOA search time;

(5)

Computing DOD and range requires $O\left({\displaystyle \sum _{i=1}^P}(N(M_{TS}-1)(2M_{TS}NK+4K^2)+K^3)\right)$. The computational complexity of this process is relatively small, so it can be ignored.

In brief, the computational complexity of the developed approach is $O\{KMNL+M^2{N}^2{K}^2(MN+K^2)+d_{\theta }\left(MK\right)!(MK-1)+1/4MNL(M+N+L)\}$.

For proving the superiority of our method, the computational complexity of the developed approach is compared with that of ESPRIT [28], Tensor-ESPRIT [30], RD-MUSIC [20], and MUSIC [27].

Table 2. Computational complexity comparison.

Method	Computational Complexity
Proposed	$O\{KMNL+M^2{N}^2{K}^2(MN+K^2)+d_{\theta }\left(MK\right)!(MK-1)$ $+1/4MNL(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\}$
Tensor-ESPRIT	$O\{2{\left(MNL\right)}^3+MNL(M+N+L)+MLK(N+K)$ $+K^3(L+M)+MN{K}^2+31K^3\}$
RD-MUSIC	$O\{4KMNL+M^2{N}^2{K}^2(MN+K^2)+d_{\theta }\left(MK\right)!(MK-1)$ $+MNL(M+N+L)\}$
MUSIC	$O\left\{L{M}^2{N}^2+4{\left(MN\right)}^3+4d_{\theta }{d}_{\phi }{d}_{r}MN(MN+1)\right\}$

shows the computational complexity comparison.

5. Simulation Results

Several numerical simulations are developed to prove the superiority of our method in parameter estimation. ESPRIT [28], tensor ESPRIT [30], and RD-MUSIC [20] are employed for performance comparison. In addition, Cramér–Rao bound (CRB) [41,42] is employed to evaluate the precision of the developed algorithm. In this part, bistatic FDA-MIMO radar with $M=18$ transmitter antennas and $N=18$ receiver antennas is considered. In the simulation, it is supposed that there are three uncorrelated targets at $({\theta }_1,{\phi }_1,r_1)=(-{50}^{\circ },-{15}^{\circ },30,000\mathrm{m})$, $({\theta }_2,{\phi }_2,r_2)=(0^{\circ },{20}^{\circ },60,000\mathrm{m})$, and $({\theta }_3,{\phi }_3,r_3)=({30}^{\circ },-{10}^{\circ },0\mathrm{m})$. To assess the precision of our method, the root mean square error (RMSE) is employed, where the RMSEs of angle and range estimation can be given by

$${\mathrm{RMSE}}_{\theta ,\phi }=\frac{1}{K}\sum _{k=1}^K\sqrt{\frac{1}{\zeta }\sum _{i=1}^{\zeta }\left\{{\left({\tilde{\theta }}_{k,i}-{\theta }_k\right)}^2+{\left({\tilde{\phi }}_{k,i}-{\phi }_k\right)}^2\right\}},$$

(55)

$${\mathrm{RMSE}}_r=\frac{1}{K}\sum _{k=1}^K\sqrt{\frac{1}{\zeta }\sum _{i=1}^{\zeta }{\left({\tilde{r}}_{k,i}-r_k\right)}^2},$$

(56)

where the results of the i-th Monte Carlo of ${\theta }_k$, ${\phi }_k$ and $r_k$ are ${\tilde{\theta }}_{k,i}$, ${\tilde{\phi }}_{k,i}$ and ${\tilde{r}}_{k,i}$. $\zeta$ is the total Monte Carlo times, $\zeta =500$.

Figure 3 and Figure 4 prove the effectiveness of the developed algorithm, where $\mathrm{SNR}=20$ dB and $L=200$. Figure 3 demonstrates the spatial spectrum of DOA. It can be seen that DOA estimation is accurately achieved through spatial spectrum searching. The spatial spectrum of the proposed algorithm is clearer and more accurate than that of MUSIC algorithm, and has more sharp peaks. We can know from Figure 4 that our method can accurately obtain the DOA, DOD, and range information of three targets. It can testify the availability of our method.

Figure 5 and Figure 6 show the comparisons of the computational complexity with the number of array elements and the number of snapshots, respectively. It is known from Figure 5 that the computational complexity of all algorithms increases with the increase in the number of array elements. The computational complexity of our algorithm is close to that of Tensor-ESPRIT and RD-MUSIC, and far lower than that of MUSIC. Figure 6 shows that the computational complexity of our method is lower than that of RD-MUSIC and slightly higher than that of Tensor-ESPRIT. In general, since our method is a real-valued operation, the computational complexity of our method is lower than that of RD-MUSIC, which is close to that of Tensor-ESPRIT algorithm.

For exploring the impact of different SNR on the precision of the algorithm, the first group of comparative tests is proposed, where $L=50$. Figure 7, Figure 8 and Figure 9 show the comparison of DOA-DOD-range estimations of different methods. From this set of comparative experiments, we can see that our method has superior precision than other algorithms in DOA, DOD, and range estimation, and nearer to the CRB. This is because the developed algorithm takes advantage of the multi-dimensional structure of data, and can obtain accurate estimations through spectral peak searching. On the one hand, the accuracy of our method in DOA and DOD estimation is very close to that of the RD-MUSIC, but the computational redundancy is reduced. On the other hand, in terms of range estimation, since the subarray division can more accurately decouple the range information and DOD information, more accurate range estimation can be obtained. Moreover, the accuracy of tensor ESPRIT is more accurate than that of ESPRIT due to the use of multi-dimensional structure characteristics.

The second group of experiments confirms the advantage of our algorithm under different snapshots, where $\mathrm{SNR}=5$ dB. From Figure 10, Figure 11 and Figure 12 we can know that DOA-DOD-range estimation accuracy of the developed algorithm are higher to the other algorithms. Moreover, the estimation accuracy of the developed algorithm is more stable in the case of low snapshots. It is because our algorithm utilizes high-dimensional data to improve the performance under small snapshots. In terms of DOA and DOD estimation, the performance of our method is slightly better than that of the RD-MUSIC. However, in range estimation, the accuracy of our method is higher to the RD-MUSIC because the subarray division solves the coupling problem of the transmitter. Moreover, the precision of tensor ESPRIT algorithm is obviously higher than ESPRIT algorithm in small snapshots. With the increase in snapshot, the precision gap between the two methods will gradually decrease. That is to say, the accuracy of the two algorithms is very close in high snapshots. Only the precision of our method is nearer to CRB curve. This demonstrates the superiority of our method.

Another criterion utilized to access the precision of the developed algorithm is probability of the successful detection (PSD), where PSD is given by

$$\mathrm{PSD}=\frac{D}{\zeta }\times 100\%,$$

(57)

where D is the number of times of obtaining the correct estimation result. If the angle error is lower than $0.1^{\circ }$ and the range error is lower than $0.1$ km, it is the correct estimation result.

For further exploring the estimation accuracy of our method, the concept of PSD is introduced to design the third group of experiments, where $L=50$. It can be seen from Figure 13, Figure 14 and Figure 15 that in DOA estimation, DOD estimation and range estimation, the PSD of the developed method is higher than the other methods at the same SNR. This means that the estimation accuracy of the developed algorithm is the highest at the same SNR. In addition, when $\mathrm{SNR}=0$ dB, PSD of our algorithm can reach $100\%$ in DOA estimation. It shows that our method can obtain more accurate results at low SNR. It is worth noting that the proposed method has less advantages than the RD-MUSIC algorithm in angle estimation, but has greater advantages in range estimation. In conclusion, the estimation accuracy of the developed approach is higher than the other methods.

6. Conclusions

In this paper, a joint DOA-DOD-range estimation algorithm with low computational complexity based on real-valued HOSVD is developed for bistatic FDA-MIMO radar. The developed algorithm can realize high-precision DOA-DOD-range estimation with low computational complexity. By employing unitary transformation technique, our method converts data into real-valued data. Then our algorithm utilizes HOSVD to estimate the signal subspace and construct a two-dimensional spatial spectral function. Then, the spatial spectral function is transformed into one-dimension by Lagrange algorithm, and the DOA estimation is obtained. The proposed approach divides the transmitting array into several subarrays, which decouples the DOD and range information. Finally, the decoupled phase is employed to achieve DOD and range estimations. Our method preserves the original multidimensional structure of the data and effectively improves the estimation precision. Numerical simulations demonstrate the preponderance of our method.