An Exact Model-Based Method for Near-Field Sources Localization with Bistatic MIMO System

Abstract

In this paper, we propose an exact model-based method for near-field sources localization with a bistatic multiple input, multiple output (MIMO) radar system, and compare it with an approximated model-based method. The aim of this paper is to propose an efficient way to use the exact model of the received signals of near-field sources in order to eliminate the systematic error introduced by the use of approximated model in most existing near-field sources localization techniques. The proposed method uses parallel factor (PARAFAC) decomposition to deal with the exact model. Thanks to the exact model, the proposed method has better precision and resolution than the compared approximated model-based method. The simulation results show the performance of the proposed method.

1. Introduction

Sources localization has been an important field of research for several decades. It is widely used in radar, underwater sources localization, acoustics, medicine, robotics, etc. The sources can be classified as near and far fields. Because of the wide range of applications, most of the research works [1,2,3,4,5] are dedicated to far-field sources localization. However, near-field sources localization has some important applications, like airport security control, ground penetration radar, phonocardiography, and many more.

Most of the existing near-field sources localization techniques [6,7,8,9,10,11,12,13,14,15,16] are based on an approximated model. In practice, a near-field point source has a spherical wavefront [6], which implies a nonlinear model. The wavefront of a near-field source is usually approximated as quadric (quadratic surface) to reduce the complexity of the model [6,9]. However, the use of this approximation results in a systematic error, which inevitably deteriorates the accuracy of the estimation. The systematic error is like an offset added to the actual source position, which increases when the target gets close to the antenna array [6].

In recent years, multiple input, multiple output (MIMO) radar has drawn a lot of attention. The advantages and limitations of MIMO radar have been well summarized in [2,3]. Based on the placement and configuration of the antennas, MIMO radar systems can broadly be classified as distributed or colocated. MIMO radar with colocated antennas can further be classified as monostatic, bistatic, and multistatic. In a bistatic MIMO radar system, the transmitting and receiving arrays are separated by a large distance, but the antennas in each array are kept close to each other (colocated) as compared to the distances between targets and the arrays. When the same array is used as transmitter and receiver, the system is monostatic. The directions of arrival and departure are different in the case of a bistatic MIMO system, but are equal for a monostatic MIMO radar system. If the distance between the transmitting and receiving arrays of a bistatic MIMO system is negligible compared to the range of targets, it can be considered as a pseudo-monostatic MIMO system. The work by Guo et al. [10] provides a subspace-based near-field sources localization method for a pseudo-monostatic MIMO radar system. A near-field sources localization method based on an approximated model with a bistatic MIMO system was proposed in [15]. Recently, in [17], a method based on the exact model of the received signal has been proposed to locate near-field targets using a bistatic MIMO system composed of linearly-aligned transmitting and receiving arrays. This paper is an improvement and an extension of the method in [17]. The major differences between the two are:

• The method in [17] is specific to the linearly-aligned transmitting and receiving arrays, whereas this paper deals with any configuration for the transmitting and receiving uniform linear arrays (ULAs) (3D situation).
• Due to the linearly-aligned transmitting and receiving arrays, the cost function in [17] has only two variables. However, the generalized 3D configuration in this paper results in a three-variable cost function which is much more difficult to deal with. Thus, in this paper, we propose a better and more efficient approach based on an overdetermined system of linear equations.

In [15], four parameters—namely, the angle of arrival, the angle of departure of a target, and the distances (ranges) from the target to the transmitting and receiving arrays—are used to localize the target, but there are some redundancies because three coordinates are sufficient to define the position of a target. Therefore, in this paper, we use Cartesian coordinates to formulate the signal model and express the localization error.

There are many existing methods to localize sources from an array of sensors, such as Capon’s method, multiple signal classification (MUSIC), estimation of signal parameter via rotational invariance techniques (ESPRIT), propagator method, and tensor decomposition method [1,5]. Among the methods listed above, the tensor decomposition method directly estimates the whole directional matrix instead of the directional parameters [18], which facilitates the estimation of the directional parameters from a nonlinear model such as the exact model in near-field situation. Consequently, the proposed exact model-based method uses the tensor decomposition. Tensor-based models and techniques are well adapted to MIMO radar because tensors allow coping with large systems (three or more dimensions). The received signal in the case of a bistatic MIMO system can be organized as a three-way tensor. Three-way tensors have attracted a lot of attention because they are the simplest form of tensor after a matrix and can be decomposed into unique factors, contrary to a matrix. Kruskal [19] provides a detailed study of the rank and uniqueness in the decomposition of a three-way tensor. The tensor decomposition has already been used for multiple far-field sources localization with bistatic MIMO radar [5].

There exist many tensor decomposition techniques, such as Tucker, parallel factor analysis (PARAFAC), and block component decomposition [20]. PARAFAC is often used in array signal processing thanks to its uniqueness in the decomposition of tensors under some mild conditions [18]. Thus, in this paper, we select PARAFAC to decompose the three-way tensor of the received signal to obtain the directional matrices of arrival and departure. From the existing work on the application of PARAFAC to the localization of targets, we can observe that it is mainly proposed for far-field target localization. In this paper, we extend it to the near-field situation. Once the directional matrices are estimated, an optimization method can be used to obtain the directional parameters.

To summarize, this paper focuses on the use of an exact model of the received signals of near-field sources to get better performance than the existing approximated model-based techniques for near-field sources localization with a bistatic MIMO system. Due to the nonlinear nature of the exact model, the PARAFAC decomposition is used, and an optimization technique is developed to efficiently solve this problem.

The remainder of the paper is organized as follows. In Section 2, a detailed signal model is constructed for a bistatic MIMO radar system based on the spherical wavefront of an incoming wave by taking the exact propagation model in the near-field situation into account. Section 3 provides a short presentation of the method proposed in [15]. In Section 4, the proposed method is described. Finally, some simulation results are presented to compare the performance of the proposed method with the method presented in [15], followed by some discussion and conclusions.

Notations

In the following, a bold lower case character (e.g., $\mathit{a}$) represents a vector, whereas a bold upper case character (e.g., $\mathit{A}$) denotes a matrix. A tensor is denoted by a bold upper case calligraphic font (e.g., $\mathcal{Y}$). ${[•]}^T$, ${[•]}^{+}$, and ${∥•∥}_F$ represent, respectively, the transpose, left pseudo-inverse, and Frobenius norm of a matrix or vector. ⊡ is the Khatri–Rao (column-wise Kronecker) product operator. The cardinal number of a set is represented by $\mathfrak{c}(•)$. $\angle (•)$ stands for the principal value of the angle (or argument) of a complex number. $\mathfrak{D}\left\{\mathit{a}\right\}$ represents the diagonal matrix with all the components of vector $\mathit{a}$ as its diagonal elements. $\mathcal{E}\{•\}$ is the expected value.

2. Signal Model

Let P be the number of narrow-band stationary point sources in the near-field region of a bistatic MIMO system with ULAs. In the following, M and N represent, respectively, the number of antennas in the transmitting and receiving arrays of the bistatic MIMO system.

For such a bistatic MIMO system, the L samples of the received matched signal in the presence of P stationary point sources can be written as [4]

$$\begin{array}{cc}\hfill {\mathit{Y}}_{\mathit{M}}& =\left({\mathit{A}}_{\mathit{e}}⊡{\mathit{A}}_{\mathit{r}}\right){\mathit{S}}^T+{\mathit{W}}_{\mathit{M}}\hfill \end{array}$$

(1)

where ${\mathit{A}}_{\mathit{e}}\in {\mathbb{C}}^{M\times P}$ and ${\mathit{A}}_{\mathit{r}}\in {\mathbb{C}}^{N\times P}$ contain the directional vectors of departure and arrival, respectively, $\mathit{S}\in {\mathbb{C}}^{L\times P}$ is the matrix of the complex-valued reflection coefficients of targets, and ${\mathit{W}}_{\mathit{M}}\in {\mathbb{C}}^{MN\times L}$ is an additive noise matrix composed of spatially- and temporally-independent elements, and each element is a zero mean Gaussian random variable with variance ${\sigma }^2$. The reflection coefficients are assumed to be different for each target and randomly changing with each sample. In other words, we consider a Swerling model II case, which makes $\mathit{S}$ a full rank matrix [21]. The pth columns of ${\mathit{A}}_{\mathit{e}}$ and ${\mathit{A}}_{\mathit{r}}$—denoted by ${\mathit{a}}_{{\mathit{e}}_p}$ and ${\mathit{a}}_{{\mathit{r}}_p}$, respectively—are given by

$$\begin{array}{c}\hfill {\mathit{a}}_{{\mathit{e}}_p}={\left[a_{e_{(1-m_o,p)}},\cdots ,1,\cdots ,a_{e_{(M-m_o,p)}}\right]}^T\end{array}$$

(2)

and

$$\begin{array}{c}\hfill {\mathit{a}}_{{\mathit{r}}_p}={\left[a_{r_{(p,1-n_o)}},\cdots ,1,\cdots ,a_{r_{(p,N-n_o)}}\right]}^T\end{array}$$

(3)

where $m_o$ and $n_o$ are the indexes of the reference elements of the transmitting and receiving arrays, respectively, $a_{e_{(m,p)}}=e^{-j2\pi {\delta }_{e_{(m,p)}}/\lambda }$, and $a_{r_{(p,n)}}=e^{-j2\pi {\delta }_{r_{(p,n)}}/\lambda }$. $m\in \{1-m_o,2-m_o,\cdots ,M-m_o\}$ and $n\in \{1-n_o,2-n_o,\cdots ,N-n_o\}$ are the relative indexes of the respective arrays. $\lambda$ is the wavelength of the carrier. ${\delta }_{e_{(m,p)}}$ is the difference between the distance traveled by the transmitted signal from the mth transmitting antenna to the pth target and the distance traveled by the transmitted signal from the 0th transmitting antenna to the pth target, which can be expressed as

$$\begin{array}{c}\hfill {\delta }_{e_{(m,p)}}=\sqrt{{\rho }_{e_p}^2+m^2{d}_e^2-2m{d}_e{\rho }_{e_p}cos\left({\theta }_{e_p}\right)}-{\rho }_{e_p}\end{array}$$

(4)

where ${\rho }_{e_p}$ and ${\theta }_{e_p}$ are respectively the range and angle of departure of the pth target with respect to the reference transmitting antenna indexed by $m_o$, and $d_e$ is the inter-element spacing in the transmitting array. Similarly, ${\delta }_{r_{(p,n)}}$ is the difference between the distance traveled by the reflected signal from the pth target to the nth receiving antenna and the distance traveled by the reflected signal from the pth target to the 0th receiving antenna, which can be expressed as

$$\begin{array}{c}\hfill {\delta }_{r_{(p,n)}}=\sqrt{{\rho }_{r_p}^2+n^2{d}_r^2-2n{d}_r{\rho }_{r_p}cos\left({\theta }_{r_p}\right)}-{\rho }_{r_p}\end{array}$$

(5)

where ${\rho }_{r_p}$ and ${\theta }_{r_p}$ are, respectively, the range and angle of arrival of the pth target with respect to the reference receiving antenna indexed by $n_o$ and $d_r$ is the inter-element spacing in the receiving array [6].

${\mathit{Y}}_{\mathit{M}}$ in Equation (1) can be considered as a block matrix, ${\mathit{Y}}_{\mathit{M}}={\left[{\breve{\mathit{Y}}}_{1-m_o}^T,{\breve{\mathit{Y}}}_{2-m_o}^T,\cdots ,{\breve{\mathit{Y}}}_{M-m_o}^T\right]}^T$. The mth sub-matrix of ${\mathit{Y}}_{\mathit{M}}$ (i.e., ${\breve{\mathit{Y}}}_m\in {\mathbb{C}}^{N\times L}$) can be expressed as

$${\breve{\mathit{Y}}}_m={\mathit{A}}_{\mathit{r}}{\mathit{D}}_m{\mathit{S}}^T+{\breve{\mathit{W}}}_m$$

(6)

where ${\mathit{D}}_m=\mathfrak{D}\left\{\left[a_{e_{(m,1)}},a_{e_{(m,2)}},\cdots ,a_{e_{(m,P)}}\right]\right\}$ and ${\breve{\mathit{W}}}_m$ is the corresponding noise sub-matrix.

3. Approximated Model-Based Method Proposed in [15]

Most of the existing near-field sources localization techniques [6,7,8,9,11,14] use an approximated model, and ULA is often used in the approximated model-based methods. The approximated path differences—which are the second-order Taylor approximations of Equations (4) and (5)—can respectively be written as [6]

$${\tilde{\delta }}_{e_{(m,p)}}=-m{d}_{e}cos\left({\theta }_{e_p}\right)+\frac{m^2{d}_e^2}{2{\rho }_{e_p}}{sin}^2\left({\theta }_{e_p}\right)$$

(7)

and

$${\tilde{\delta }}_{r_{(p,n)}}=-n{d}_{r}cos\left({\theta }_{r_p}\right)+\frac{n^2{d}_r^2}{2{\rho }_{r_p}}{sin}^2\left({\theta }_{r_p}\right).$$

(8)

In [15], a subspace-based method is used to estimate four parameters: two ranges (${\rho }_{e_p}$ and ${\rho }_{r_p}$) and two directional angles (${\theta }_{e_p}$ and ${\theta }_{r_p}$) of a near-field target by using a bistatic MIMO system with inter-element spacing of $\lambda /4$ in each ULA. $M\in \{2\mu +1:\mu \in \mathbb{N},\mu >1\}$, $N\in \{2\breve{\mu }+1:\breve{\mu }\in \mathbb{N}\}$, $m_o=(M+1)/2$, $n_o=(N+1)/2$, $d_e\le \lambda /4$, and $d_r\le \lambda /4$ are the necessary conditions of [15].

In an approximated model-based method like [15], ${\mathit{A}}_{\mathit{e}}$ and ${\mathit{A}}_{\mathit{r}}$ are assumed to be constructed by ${\tilde{\delta }}_{e_{(m,p)}}$ and ${\tilde{\delta }}_{r_{(p,n)}}$, respectively. Therefore, in this case, ${\mathit{D}}_m$ in Equation (6) can be expressed as

$${\mathit{D}}_m=\mathfrak{D}\left\{\left[e^{j\left(m{\omega }_{e_1}-m^2{\varphi }_{e_1}\right)},e^{j\left(m{\omega }_{e_2}-m^2{\varphi }_{e_2}\right)},\cdots ,e^{j\left(m{\omega }_{e_P}-m^2{\varphi }_{e_P}\right)}\right]\right\}$$

(9)

where ${\omega }_{e_p}=2\pi d_{e}cos\left({\theta }_{e_p}\right)/\lambda$ and ${\varphi }_{e_p}=\pi d_e^2{sin}^2\left({\theta }_{e_p}\right)/\left(\lambda {\rho }_{e_p}\right)$. In [15], four cross-covariance matrices between ${\breve{\mathit{Y}}}_m$ for $m\in \{-2,-1,1,2\}$ and ${\breve{\mathit{Y}}}_0$ are constructed. The eigenvalues of ${\mathit{R}}_{-2}{\mathit{R}}_{-1}^{+}$ and ${\mathit{R}}_2{\mathit{R}}_1^{+}$ are used to get ${\rho }_{e_p}$ and ${\theta }_{e_p}$, and their eigenvectors are used to obtain ${\rho }_{r_p}$ and ${\theta }_{r_p}$, where ${\mathit{R}}_m=\mathcal{E}\{{\breve{\mathit{Y}}}_m{\breve{\mathit{Y}}}_0^H\}$. More details can be found in [15].

4. Proposed Exact Model-Based Position Estimation Method

Every element of ${\mathit{Y}}_{\mathit{M}}$ in Equation (1) is associated with three parameters related, respectively, to the receiving antenna, transmitting antenna, and time sample. Therefore, ${\mathit{Y}}_{\mathit{M}}$ can be rearranged like a three-way tensor $\mathcal{Y}\in {\mathbb{C}}^{N\times M\times L}$, as shown in Figure 1. Creating a tensor out of lower dimensional data is known as tensorization [20].

PARAFAC decomposition of tensor $\mathcal{Y}$ is used to get the estimates of ${\mathit{A}}_{\mathit{r}}$, ${\mathit{A}}_{\mathit{e}}$, and $\mathit{S}$ matrices [5]. Tensor operations are usually performed in its equivalent matrix form [5,22,23]. The process of creating a matrix out of a tensor is known as matricization [20]. Like ${\mathit{Y}}_{\mathit{M}}$, $\mathcal{Y}$ can be matricized into the following two additional matrices

$${\mathit{Y}}_{\mathit{L}}=\left(\mathit{S}⊡{\mathit{A}}_{\mathit{e}}\right){{\mathit{A}}_{\mathit{r}}}^T+{\mathit{W}}_{\mathit{L}}$$

(10)

and

$${\mathit{Y}}_{\mathit{N}}=\left({\mathit{A}}_{\mathit{r}}⊡\mathit{S}\right){{\mathit{A}}_{\mathit{e}}}^T+{\mathit{W}}_{\mathit{N}}.$$

(11)

According to the least squares principle, the following objective functions can be written from Equations (1), (10) and (11)

$$\widehat{\mathit{S}}=\underset{\mathit{S}}{arg\; min}\left\{{∥{\mathit{Y}}_{\mathit{M}}-\left({\mathit{A}}_{\mathit{e}}⊡{\mathit{A}}_{\mathit{r}}\right){\mathit{S}}^T∥}_F^2\right\},$$

(12)

$${\widehat{\mathit{A}}}_{\mathit{r}}=\underset{{\mathit{A}}_{\mathit{r}}}{arg\; min}\left\{{∥{\mathit{Y}}_{\mathit{L}}-\left(\mathit{S}⊡{\mathit{A}}_{\mathit{e}}\right){{\mathit{A}}_{\mathit{r}}}^T∥}_F^2\right\},$$

(13)

and

$${\widehat{\mathit{A}}}_{\mathit{e}}=\underset{{\mathit{A}}_{\mathit{e}}}{arg\; min}\left\{{∥{\mathit{Y}}_{\mathit{N}}-\left({\mathit{A}}_{\mathit{r}}⊡\mathit{S}\right){{\mathit{A}}_{\mathit{e}}}^T∥}_F^2\right\}$$

(14)

where ${\widehat{\mathit{A}}}_{\mathit{r}}$, ${\widehat{\mathit{A}}}_{\mathit{e}}$, and $\widehat{\mathit{S}}$ denote the estimated values of ${\mathit{A}}_{\mathit{r}}$, ${\mathit{A}}_{\mathit{e}}$, and $\mathit{S}$ respectively.

The trilinear alternating least squares (TALS) algorithm is a classical method to minimize the above objective functions [5,22,23]. Least squares estimates of Equations (12)–(14) are given by

$$\widehat{\mathit{S}}={\left[{\left({\mathit{A}}_{\mathit{e}}⊡{\mathit{A}}_{\mathit{r}}\right)}^{+}{\mathit{Y}}_{\mathit{M}}\right]}^T,$$

(15)

$${\widehat{\mathit{A}}}_{\mathit{r}}={\left[{\left(\mathit{S}⊡{\mathit{A}}_{\mathit{e}}\right)}^{+}{\mathit{Y}}_{\mathit{L}}\right]}^T,$$

(16)

and

$${\widehat{\mathit{A}}}_{\mathit{e}}={\left[{\left({\mathit{A}}_{\mathit{r}}⊡\mathit{S}\right)}^{+}{\mathit{Y}}_{\mathit{N}}\right]}^T.$$

(17)

In the TALS algorithm, Equations (15)–(17) are alternatively updated with the new values of ${\widehat{\mathit{A}}}_{\mathit{r}}$, ${\widehat{\mathit{A}}}_{\mathit{e}}$, and $\widehat{\mathit{S}}$ until a stopping criteria is met. $∥{\mathit{Y}}_{\mathit{M}}-{\left({\mathit{A}}_{\mathit{e}}⊡{\mathit{A}}_{\mathit{r}}\right){\mathit{S}}^T∥}_F^2<{ϵ}_{\mathrm{tol}}$ is often used as the stopping condition, where ${ϵ}_{\mathrm{tol}}$ is the tolerance. In practice, the algorithm given in [24] is used for PARAFAC decomposition, which uses compression, line search, normalization, etc. to accelerate the TALS method.

According to [19], the matrices obtained by PARAFAC decomposition of a three-way tensor are scaled and permuted. The permutation has no impact because the matrices’ columns are paired. However, in the proposed method, the scaling factor must be removed by dividing all the elements of the directional vectors with their corresponding reference elements.

To define the Cartesian coordinates of a target, let us assume a general configuration of bistatic MIMO system as shown in Figure 2. In the case of a ULA, the unit vector along the array and the position vector of the reference antenna of the corresponding array are sufficient to obtain the position vectors of the remaining antennas of that array. In the figure, ${\mathit{e}}_o$ and ${\mathit{r}}_o$ are the position vectors of the reference transmitting and receiving antennas, respectively, with respect to the origin of the Cartesian coordinate system, and ${\mathit{d}}_{{\mathit{c}}_{\mathit{e}}}$ and ${\mathit{d}}_{{\mathit{c}}_{\mathit{r}}}$ are the unit vectors along the transmitting and receiving arrays, respectively. ${\mathit{t}}_p={[x_{t_p},y_{t_p},z_{t_p}]}^T$ represents the position vector of the pth target. In 3D space, the range and directional angle of a target with respect to a linear array make a circle related to the base of a cone with the range as its slant height and the directional angle as its half angle. In the bistatic case, we have two such circles (shown in Figure 2). The target is located at the intersection of these two circles. In the figure, ${\mathit{\nu }}_{{\mathit{e}}_p}$ and ${\mathit{\nu }}_{{\mathit{r}}_p}$ are unit vectors on the planes of the respective circles. The parametric equations of the circles can be written as

$${\mathit{\psi }}_{{\mathit{e}}_p}\left({\phi }_e\right)={\rho }_{e_p}sin\left({\theta }_{e_p}\right)\left[cos\left({\phi }_e\right){\mathit{\nu }}_{{\mathit{e}}_p}+sin\left({\phi }_e\right){\mathit{d}}_{{\mathit{c}}_{\mathit{e}}}\times {\mathit{\nu }}_{{\mathit{e}}_p}\right]+{\rho }_{e_p}cos\left({\theta }_{e_p}\right){\mathit{d}}_{{\mathit{c}}_{\mathit{e}}}+{\mathit{e}}_o$$

(18)

and

$${\mathit{\psi }}_{{\mathit{r}}_p}\left({\phi }_r\right)={\rho }_{r_p}sin\left({\theta }_{r_p}\right)\left[cos\left({\phi }_r\right){\mathit{\nu }}_{{\mathit{r}}_p}+sin\left({\phi }_r\right){\mathit{d}}_{{\mathit{c}}_{\mathit{r}}}\times {\mathit{\nu }}_{{\mathit{r}}_p}\right]+{\rho }_{r_p}cos\left({\theta }_{r_p}\right){\mathit{d}}_{{\mathit{c}}_{\mathit{r}}}+{\mathit{r}}_o$$

(19)

where × denotes the cross-product operation between two vectors; ${\mathit{\psi }}_{{\mathit{e}}_p}\left({\phi }_e\right)$ and ${\mathit{\psi }}_{{\mathit{r}}_p}\left({\phi }_r\right)$ are the position vectors of a point on the respective circles at ${\phi }_e$ and ${\phi }_r$, respectively. The equation parameters ${\phi }_e$ and ${\phi }_r$ independently vary from 0 to $2\pi$ rad to completely sweep the respective circles.

The ranges and directional angles can be expressed in terms of the Cartesian coordinates as ${\rho }_{e_p}={∥{\mathit{t}}_p-{\mathit{e}}_o∥}_F$, ${\rho }_{r_p}={∥{\mathit{t}}_p-{\mathit{r}}_o∥}_F$, ${\theta }_{e_p}=arccos\left[{\left({\mathit{t}}_p-{\mathit{e}}_o\right)}^T{\mathit{d}}_{{\mathit{c}}_{\mathit{e}}}/{\rho }_{e_p}\right]$, and ${\theta }_{r_p}=arccos\left[{\left({\mathit{t}}_p-{\mathit{r}}_o\right)}^T{\mathit{d}}_{{\mathit{c}}_{\mathit{r}}}/{\rho }_{r_p}\right]$. Thus, according to Equations (2)–(5), ${\mathit{a}}_{{\mathit{e}}_p}$ and ${\mathit{a}}_{{\mathit{r}}_p}$ can respectively be determined by ${\mathit{t}}_p$ as

$$a_{e_{(m,p)}}=e^{-j2\pi \left(\sqrt{{∥{\mathit{t}}_p-{\mathit{e}}_o∥}_F^2+m^2{d}_e^2-2m{d}_e{\left({\mathit{t}}_p-{\mathit{e}}_o\right)}^T{\mathit{d}}_{{\mathit{c}}_{\mathit{e}}}}-{∥{\mathit{t}}_p-{\mathit{e}}_o∥}_F\right)/\lambda }$$

(20)

and

$$a_{r_{(p,n)}}=e^{-j2\pi \left(\sqrt{{∥{\mathit{t}}_p-{\mathit{r}}_o∥}_F^2+n^2{d}_r^2-2n{d}_r{\left({\mathit{t}}_p-{\mathit{r}}_o\right)}^T{\mathit{d}}_{{\mathit{c}}_{\mathit{r}}}}-{∥{\mathit{t}}_p-{\mathit{r}}_o∥}_F\right)/\lambda }.$$

(21)

Then, a direct approach to estimate ${\mathit{t}}_p$ could be the minimization of the following cost function

$$\begin{array}{c}\hfill {\widehat{\mathit{t}}}_p=\underset{{\mathit{t}}_p}{arg\; min}\left\{{∥{\widehat{\mathit{a}}}_{{\mathit{e}}_p}/{\widehat{a}}_{e_{(0,p)}}-{\mathit{a}}_{{\mathit{e}}_p}∥}_F^2+{∥{\widehat{\mathit{a}}}_{{\mathit{r}}_p}/{\widehat{a}}_{r_{(p,0)}}-{\mathit{a}}_{{\mathit{r}}_p}∥}_F^2\right\}\end{array}$$

(22)

where ${\widehat{\mathit{a}}}_{{\mathit{e}}_p}$ and ${\widehat{\mathit{a}}}_{{\mathit{r}}_p}$ are the estimated directional vectors obtained from the PARAFAC decomposition, and ${\widehat{a}}_{e_{(0,p)}}$ and ${\widehat{a}}_{r_{(p,0)}}$ are their respective reference elements used here to remove the scaling factor in the decomposition of the received signal tensor.

Even though a near-field region occupies a finite space, minimizing Equation (22) by using grid search or Newton’s method is computationally expensive. Therefore, we choose an indirect method in which we estimate the ranges and directional angles, followed by the estimation of the coordinates.

Rearranging Equations (4) and (5), we can obtain

$$\begin{array}{cc}\hfill 2m{d}_e{\rho }_{e_p}cos\left({\theta }_{e_p}\right)+2{\widehat{\delta }}_{e_{(m,p)}}{\rho }_{e_p}& =m^2{d}_e^2-{\widehat{\delta }}_{e_{(m,p)}}^2\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill 2n{d}_r{\rho }_{r_p}cos\left({\theta }_{r_p}\right)+2{\widehat{\delta }}_{r_{(p,n)}}{\rho }_{r_p}& =n^2{d}_r^2-{\widehat{\delta }}_{r_{(p,n)}}^2\hfill \end{array}$$

(24)

where ${\widehat{\delta }}_{e_{(m,p)}}$ and ${\widehat{\delta }}_{r_{(p,n)}}$ are the estimated path differences which can be directly obtained from the estimated directional vectors as follows

$${\widehat{\delta }}_{e_{(m,p)}}=-\lambda \left[\mathfrak{U}\left\{\angle \left({\widehat{a}}_{e_{(m,p)}}\right)\right\}\right.-\left.\mathfrak{U}\left\{\angle \left({\widehat{a}}_{e_{(0,p)}}\right)\right\}\right]/2\pi$$

(25)

and

$${\widehat{\delta }}_{r_{(p,n)}}=-\lambda \left[\mathfrak{U}\left\{\angle \left({\widehat{a}}_{r_{(p,n)}}\right)\right\}\right.-\left.\mathfrak{U}\left\{\angle \left({\widehat{a}}_{r_{(p,0)}}\right)\right\}\right]/2\pi$$

(26)

where $\mathfrak{U}\left\{•\right\}$ represents the unwrapped value of the argument [25]. Equations (25) and (26) can be described as follows:

• Get the directional vectors ${\widehat{\mathit{a}}}_{{\mathit{e}}_p}$ and ${\widehat{\mathit{a}}}_{{\mathit{r}}_p}$ from the PARAFAC decomposition.
• Extract the arguments of all the components of ${\widehat{\mathit{a}}}_{{\mathit{e}}_p}$ and ${\widehat{\mathit{a}}}_{{\mathit{r}}_p}$.
• Unwrap the phase vectors obtained from Step 2.
• Subtract the unwrapped phase corresponding to ${\widehat{a}}_{e_{(0,p)}}$ and ${\widehat{a}}_{r_{(p,0)}}$ from all the components of the unwrapped phase vector of ${\widehat{\mathit{a}}}_{{\mathit{e}}_p}$ and ${\widehat{\mathit{a}}}_{{\mathit{r}}_p}$, respectively.
• Divide each component of the normalized phase vectors obtained from the above step by $-2\pi /\lambda$ to get ${\widehat{\delta }}_{e_{(m,p)}}$ and ${\widehat{\delta }}_{r_{(p,n)}}$.

In practice, $M\ge 2$ and $N\ge 2$; therefore, Equations (23) and (24) can be considered as an overdetermined system of linear equations in ${\rho }_{e_p}cos\left({\theta }_{e_p}\right)$ and ${\rho }_{e_p}$ and ${\rho }_{r_p}cos\left({\theta }_{r_p}\right)$ and ${\rho }_{r_p}$, respectively, which can be solved by the total least squares method [26]. Let ${[u_{1_p},u_{2_p},u_{3_p}]}^T$ be the right-singular-vector associated with the smallest singular value of the following matrix formed by the coefficients of Equation (23)

$$\left[\begin{array}{ccc}2(1-m_o)d_e& 2{\widehat{\delta }}_{e_{(1-m_o,p)}}& {(1-m_o)}^2{d}_e^2-{\widehat{\delta }}_{e_{(1-m_o,p)}}^2\\ 2(2-m_o)d_e& 2{\widehat{\delta }}_{e_{(2-m_o,p)}}& {(2-m_o)}^2{d}_e^2-{\widehat{\delta }}_{e_{(2-m_o,p)}}^2\\ ⋮& ⋮& ⋮\\ 2(M-m_o)d_e& 2{\widehat{\delta }}_{e_{(M-m_o,p)}}& {(M-m_o)}^2{d}_e^2-{\widehat{\delta }}_{e_{(M-m_o,p)}}^2\end{array}\right].$$

(27)

The estimated range and angle of departure can respectively be computed by ${\widehat{\rho }}_{e_p}=-u_{2_p}/u_{3_p}$ and ${\widehat{\theta }}_{e_p}=arccos\left(u_{1_p}/u_{2_p}\right)$. Similarly, let ${[v_{1_p},v_{2_p},v_{3_p}]}^T$ be the right-singular-vector associated with the smallest singular value of the following matrix formed by the coefficients of Equation (24)

$$\left[\begin{array}{ccc}2(1-n_o)d_r& 2{\widehat{\delta }}_{r_{(p,1-n_o)}}& {(1-n_o)}^2{d}_r^2-{\widehat{\delta }}_{r_{(p,1-n_o)}}^2\\ 2(2-n_o)d_r& 2{\widehat{\delta }}_{r_{(p,2-n_o)}}& {(2-n_o)}^2{d}_r^2-{\widehat{\delta }}_{r_{(p,2-n_o)}}^2\\ ⋮& ⋮& ⋮\\ 2(N-n_o)d_r& 2{\widehat{\delta }}_{r_{(p,N-n_o)}}& {(N-n_o)}^2{d}_r^2-{\widehat{\delta }}_{r_{(p,N-n_o)}}^2\end{array}\right].$$

(28)

The estimated range and angle of arrival can respectively be computed by ${\widehat{\rho }}_{r_p}=-v_{2_p}/v_{3_p}$ and ${\widehat{\theta }}_{r_p}=arccos\left(v_{1_p}/v_{2_p}\right)$.

The estimated ranges and directional angles can be used in Equations (18) and (19) to construct the parametric equations of the circles. As mentioned before, the required coordinates are at the intersection of these circles. However, due to the estimation error and noise, the circles may not intersect; thus, the following minimization problem can be solved:

$$\begin{array}{c}\hfill ({\widehat{\phi }}_e,{\widehat{\phi }}_r)=\underset{({\phi }_e,{\phi }_e)}{arg\; min}\left\{{∥{\mathit{\psi }}_{{\mathit{e}}_p}\left({\phi }_e\right)-{\mathit{\psi }}_{{\mathit{r}}_p}\left({\phi }_r\right)∥}_F^2\right\}.\end{array}$$

(29)

A coarse solution of Equation (29) can be calculated by exhaustive grid search, and then it can be finely tuned by Newton’s method. Solving Equation (29) is less complex than solving Equation (22). Finally, the position vector of the pth target can be computed as ${\widehat{\mathit{t}}}_p={\mathit{\psi }}_{{\mathit{e}}_p}\left({\widehat{\phi }}_e\right)$, ${\widehat{\mathit{t}}}_p={\mathit{\psi }}_{{\mathit{r}}_p}\left({\widehat{\phi }}_r\right)$, or the average of these two position vectors.

Algorithm 1 provides a summary of the proposed method.

Algorithm 1 Algorithm of the proposed method.

Construct the three-way tensor $\mathcal{Y}$ from the received data. Estimate ${\mathit{A}}_{\mathit{e}}$ and ${\mathit{A}}_{\mathit{r}}$ from $\mathcal{Y}$ using PARAFAC decomposition. Use Equations (25) and (26) to obtain ${\widehat{\delta }}_{e_{(m,p)}}$ and ${\widehat{\delta }}_{r_{(p,n)}}$ from the estimated ${\mathit{A}}_{\mathit{e}}$ and ${\mathit{A}}_{\mathit{r}}$, respectively. Create the system of linear equations by substituting ${\widehat{\delta }}_{e_{(m,p)}}$ and ${\widehat{\delta }}_{r_{(p,n)}}$ in Equations (23) and (24), respectively, for all the values of m and n for each target. Separately solve each system of linear equations created in step 4 using the total least squares technique to obtain ${\widehat{\rho }}_{e_p}$, ${\widehat{\theta }}_{e_p}$, ${\widehat{\rho }}_{r_p}$, and ${\widehat{\theta }}_{r_p}$. Substitute the four estimated location parameters in Equations (18) and (19) to obtain the parametric equations of the circles, and minimize (29) to estimate ${\widehat{\phi }}_e$ and ${\widehat{\phi }}_r$. Finally, substitute ${\widehat{\phi }}_e$ and ${\widehat{\phi }}_r$ along with ${\widehat{\rho }}_{e_p}$, ${\widehat{\theta }}_{e_p}$, ${\widehat{\rho }}_{r_p}$, and ${\widehat{\theta }}_{r_p}$ in Equations (18) and (19) to get the estimated coordinates ${\widehat{\mathit{t}}}_p$ of the pth target.

5. Simulation Results

In the following simulations, the performance of the proposed method and the method in [15] is compared, with $M=5$, $N=9$, $m_o=3$, $n_o=5$, and $d_e=d_r=\lambda /4$ to satisfy the necessary requirements of [15]. Throughout the simulation, $\lambda$ is used as the unit of length. The remaining MIMO system configuration parameters are ${\mathit{e}}_o={\left[0.5\lambda ,\lambda ,1.5\lambda \right]}^T$, ${\mathit{r}}_o={\left[0.4\lambda ,-0.3\lambda ,-0.2\lambda \right]}^T$, ${\mathit{d}}_{{\mathit{c}}_{\mathit{e}}}={\left[0.3420,0.5000,-0.7956\right]}^T$, and ${\mathit{d}}_{{\mathit{c}}_{\mathit{r}}}={\left[0.8660,-0.1736,0.4689\right]}^T$, which are chosen randomly such that there exists a significant near-field region shared by both ULAs.

According to the three estimated coordinates, the root mean square error (RMSE) associated with the position estimation of the pth target is calculated as follows:

$$\begin{array}{c}\hfill {\mathrm{RMSE}}_p=\sqrt{\frac{1}{K}{\displaystyle \sum _{k=1}^K}{∥{\widehat{\mathit{t}}}_p\left(k\right)-{\mathit{t}}_p∥}_F^2}\end{array}$$

(30)

where K is the number of Monte Carlo iterations, ${\widehat{\mathit{t}}}_p\left(k\right)$ represents the estimated position at the kth iteration, and ${\mathit{t}}_p$ is the true position of the pth target.

In Figure 3, we have compared the performance of the proposed method with the method proposed in [15] with two targets at ${\left[\lambda ,\lambda ,\lambda \right]}^T$ and ${\left[2\lambda ,1.75\lambda ,1.5\lambda \right]}^T$ in the Fresnel region. The cost Function (29) has also been applied to [15] to obtain the coordinates. In addition, we have also drawn the Cramér–Rao lower bound (CRLB) in Figure 3, which can be obtained from the existing works [4,13,27]. To keep mathematical analogy with RMSE, we combine the CRLB of the three coordinates of the pth target as

$$\begin{array}{c}\hfill {CRLB}_p=\sqrt{{\sigma }_{x_p}^2+{\sigma }_{y_p}^2+{\sigma }_{z_p}^2}\end{array}$$

(31)

where ${\sigma }_{x_p}^2$, ${\sigma }_{y_p}^2$, and ${\sigma }_{z_p}^2$ denote the CRLB of the corresponding coordinates belonging to the pth target.

Figure 3 shows that the proposed method has higher precision and much better performance in terms of RMSE than that of the subspace and approximated model-based method in [15]. This performance gain comes principally from the use of the exact near-field received signal model and PARAFAC decomposition in the proposed technique.

The resolution capability of a method can be evaluated by the probability of the successful detection $P\left(\xi \right)$ of two closely-placed targets, which can be calculated as [12]:

$$\begin{array}{c}\hfill P\left(\xi \right)=\mathfrak{c}\left(\left\{k:{∥{\widehat{\mathit{t}}}_1\left(k\right)-{\mathit{t}}_1∥}_F<\xi \mathrm{and}{∥{\widehat{\mathit{t}}}_2\left(k\right)-{\mathit{t}}_2∥}_F<\xi \right\}\right)/K\end{array}$$

(32)

where $k\in \{1,2,\cdots ,K\}$ and $\xi ={∥{\mathit{t}}_1-{\mathit{t}}_2∥}_F/2$ is the half of the distance between the two targets. Figure 4 gives the probability of successful detection of two targets at different distances between the targets, which shows that the proposed method has a much better resolution power than its counterpart, even at the high signal-to-noise ratio (SNR) of 10 dB.

6. Discussion

In Figure 3, we can observe a significant gap between the RMSE corresponding to the proposed method and the method in [15]. The gain in performance of the proposed method comes from the use of the PARAFAC decomposition and the exact model of the received near-field signals. At high SNR, the method in [15] experiences a floor effect in terms of achievable RMSE performance, which clearly shows the systematic bias introduced by the approximated model. This systematic error is not discernible at low SNR because the major contribution to the estimation error comes from the noise. This bias also explains the low successful detection probability of [15], as shown in Figure 4. Because of the approximation, the location estimated by an approximated model-based method is shifted from the true location, which makes $P\left(\xi \right)$ small.

7. Conclusions

In this paper, we propose a novel technique for near-field sources localization with a bistatic MIMO system. The principal originalities of this work are the use of the exact model and PARAFAC decomposition for near-field sources localization. Thanks to the exact model of near-field sources, the proposed method has high precision and resolution. The performance of the proposed method greatly surpasses the high-resolution subspace-based method proposed in [15], which proves the importance of an exact model-based method. The proposed method also has some additional advantages with respect to the compared approximated model-based method: it works for the inter-element spacing of $\lambda /2$ without any ambiguity, M and N are not required to be odd, and any antenna can be used as the reference point.