Geometric Algebra-Based ESPRIT Algorithm for DOA Estimation

Abstract

Direction-of-arrival (DOA) estimation plays an important role in array signal processing, and the Estimating Signal Parameter via Rotational Invariance Techniques (ESPRIT) algorithm is one of the typical super resolution algorithms for direction finding in an electromagnetic vector-sensor (EMVS) array; however, existing ESPRIT algorithms treat the output of the EMVS array either as a “long vector”, which will inevitably lead to loss of the orthogonality of the signal components, or a quaternion matrix, which may result in some missing information. In this paper, we propose a novel ESPRIT algorithm based on Geometric Algebra (GA-ESPRIT) to estimate 2D-DOA with double parallel uniform linear arrays. The algorithm combines GA with the principle of ESPRIT, which models the multi-dimensional signals in a holistic way, and then the direction angles can be calculated by different GA matrix operations to keep the correlations among multiple components of the EMVS. Experimental results demonstrate that the proposed GA-ESPRIT algorithm is robust to model errors and achieves less time complexity and smaller memory requirements.

1. Introduction

Direction-of-arrival (DOA) estimation of electromagnetic (EM) signals has attracted wide attention in many communication fields, such as radar [1,2], mobile networks [3] and sonar [4]. It is clear that DOA estimation is the basic and essential part in an array signal processing system. For example, a corresponding transmitting or receiving beamformer can be designed to extract signals in the direction of interest and suppress uninteresting interference signals. The electromagnetic vector sensor (EMVS) can catch polarization-related information compared to a conventional scalar sensor, which can further improve the target resolution, anti-interference ability and detection stability for DOA estimation [5,6,7]; therefore, the research for EMVS array direction finding has become a hotspot.

With the appearance of the Long-Vector MODEL (LV-MODEL) [5] (built for EMVS), multiple researchers have proposed various DOA estimators. The existing estimators can be summarized into three categories: (1) research on DOA estimators transplanting from scalar sensor; (2) research based on special array arrangement; (3) research based on advanced mathematical tools.

In terms of transplantation, the classic subspace-based super-resolution algorithm [8] (Multiple Signal Classification—MUSIC) was transplanted to the EMVS [9,10,11] array, but the algorithms often suffer high computational complexity because of the four-dimensional parameter search for two direction angles and two additional polarization angles; therefore, Weiss [12] used the polynomial root to reduce the computational complexity to a certain extent. In addition, another subspace-based super-resolution algorithm [13,14] (Estimation of Signal Parameters via Rotational Invariance Techniques—ESPRIT) was also transplanted into the EMVS array, and realized closed-form estimation of DOA. In [15,16], authors showed that the statistical performance of the maximum likelihood and subspace-fitting algorithms based on the EMVS array are better than both MUSIC and ESPRIT, but the high calculation limits its application in actual engineering.

There are few studies based on the special array arrangement because most EMVS arrays are co-centered, leading to the mutual coupling interference and spatial information loss. In [17], a double-parallel-line EMVS array whose six components are all spatially separated achieved mutual coupling reduction to refine the DOA-finding accuracy by orders of magnitude. A triangular array [18] combined with a vector cross product and interferometric angle measurement, aimed to overcome the drawback that [17] cannot achieve two-dimensional aperture expansion. In addition, a spatial expansion method of a triangle structure [19] was proposed to provide higher-precision DOA estimation.

The traditional model for EMVS is just a linear combination of each component, which somehow locally destroy the orthogonality of the signal components [20]. Meanwhile, the heavy computational efforts and memory requirements during data processing for the DOA estimation cannot be ignored [21]. Recently, the hypercomplex has been widely studied and applied in multi-dimensional parameter estimation. Miron et al. [22] first proposed a new Quaternion Model (Q-MODEL) for the two-component EMVS array. Then, many models and algorithms based on quaternion have been proposed [23,24,25]; however, the Q-MODEL had to discard some of the original information because the quaternion only has three imaginary parts. Further, the research has extended to bi-quaternion [26] and quad-quaternion [27,28]. These quaternion-based algorithms showed higher estimation accuracy and less complexity; however, Jiang et al. [21] found that the physical interpretations of the presented quaternion-like models have not been discussed. In order to solve the problem, they derived G-MODEL [21] by Geometric Algebra (GA) formulations of Maxwell equations. The computing technology of G-MODEL not only minimizes the memory requirements and computational complexity, but also removes the correlation of noise on different antennas.

It is easy to find that the current studies utilizing hypercomplex algebra are mainly focused on the MUSIC algorithm [22,26,27,28]. In fact, MUSIC greatly suffers from a heavy computational burden for its spectrum search, while the computation of ESPRIT algorithm is cheaper, and it can automatically decouple [29]; therefore, the research in this paper extends the ESPRIT algorithm using a new mathematical tool—GA. Through the new calculation rules, the physical nature of EMVS is matched with the signal processing technology, which avoids correlation loss between different components in the previous algorithms. The major contributions of this paper are as follows.

• We incorporate the multi-dimensional consistency of GA into ESPRIT, and propose a Geometric Algebra-based ESPRIT algorithm (GA-ESPRIT) for 2D-DOA estimation.
• We use the new calculation rules of the high-dimensional algebra system to preserve the correlation among multiple components of EMVS.
• Experimental results demonstrate that the proposed GA-ESPRIT algorithm can achieve more accurate, stable and lighter DOA estimation.

The rest of this paper is organized as follows. Section 2 introduces the basics of GA and the EMVS model for narrow-band signals based on GA. Section 3 describes the proposed GA-ESPRIT in detail. Experimental results and analysis are provided in Section 4, followed by concluding remarks in Section 5.

2. Preliminaries

2.1. Fundamental of Geometric Algebra

The concept of GA [30] was proposed by David Hestenes in the 1960s, who combined Clifford Algebra with a physical geometric structure. After decades of research, GA has shown its absolute superiority in electromagnetism [31], cosmology [32], multi-channel image [33,34,35] and other physical sciences.

2.1.1. Geometric Product

The crucial product operation in GA theory is the geometric product [30]. For vectors $\mathbf{a}$ and $\mathbf{b}$, the geometric product is denoted by

$$\mathbf{a}\mathbf{b}=\mathbf{a}·\mathbf{b}+\mathbf{a}\wedge \mathbf{b},$$

(1)

where $\{·\}$ and $\{\wedge \}$ denote the inner product and the outer product, respectively.

2.1.2. Multi-Vector

Let ${\mathbb{G}}_n=C{\ell }_{n,0}$, which is the real GA of the quadratic pair $(V,Q)$ where $V={\mathbb{R}}^n$ and Q is the quadratic form of signature $(n,0)$. There is an orthogonal basis $\left\{{\mathbf{e}}_1,{\mathbf{e}}_2,\dots ,{\mathbf{e}}_n\right\}$ in ${\mathbb{R}}^n$, which generates $2^n$ basis elements of ${\mathbb{G}}_n$ via the geometric product as shown in (2):

$$\underset{k=0}{\underbrace{\left\{1\right\}}},\underset{k=1}{\underbrace{\left\{{\mathbf{e}}_i\right\}}},\underset{k=2}{\underbrace{\left\{{\mathbf{e}}_{ij},i<j\right\}}},\dots ,\underset{k=n}{\underbrace{\left\{{\mathbf{e}}_1{\mathbf{e}}_2\cdots {\mathbf{e}}_n\right\}}}$$

(2)

for $i,j=1,2,\dots ,n$.

The multi-vector A of ${\mathbb{G}}_n$ is defined as

$$\begin{array}{cc}\hfill A& =E_0\left(A\right)+\sum _{1\le i\le n}{E}_i\left(A\right){\mathbf{e}}_i+\sum _{1\le i<j\le n}{E}_{ij}\left(A\right){\mathbf{e}}_{ij}+\dots +E_{1\dots n}\left(A\right){\mathbf{e}}_{1\dots n}\hfill \\ \hfill & ={\langle A\rangle }_0+{\langle A\rangle }_1+{\langle A\rangle }_2+\dots +{\langle A\rangle }_n,\hfill \end{array}$$

(3)

where $E_i\left(A\right),E_{ij}\left(A\right),\dots ,E_{1\dots n}\left(A\right)\in \mathbb{R}$, and ${\langle A\rangle }_k$ denotes the component of A of grade k.

The reverse of multi-vector A is defined as

$$\tilde{A}=\sum _{k=0}^n{(-1)}^{k(k-1)/2}{\langle A\rangle }_k.$$

(4)

2.2. The Geometric Algebra of Euclidean 3-Space

According to the structural characteristics of EMVS, ${\mathbb{G}}_3$ is chosen to model and process the received signals [21]. The multiplication rule can be found in

Table 1. The multiplication rule in ${\mathbb{G}}_3$.

	1	${\mathbf{e}}_1$	${\mathbf{e}}_2$	${\mathbf{e}}_3$	${\mathbf{e}}_{12}$	${\mathbf{e}}_{23}$	${\mathbf{e}}_{13}$	${\mathbf{e}}_{123}$
1	1	${\mathbf{e}}_1$	${\mathbf{e}}_2$	${\mathbf{e}}_3$	${\mathbf{e}}_{12}$	${\mathbf{e}}_{23}$	${\mathbf{e}}_{13}$	${\mathbf{e}}_{123}$
${\mathbf{e}}_1$	${\mathbf{e}}_1$	1	${\mathbf{e}}_{12}$	${\mathbf{e}}_{13}$	${\mathbf{e}}_2$	${\mathbf{e}}_{123}$	${\mathbf{e}}_3$	${\mathbf{e}}_{23}$
${\mathbf{e}}_2$	${\mathbf{e}}_2$	$-{\mathbf{e}}_{12}$	1	${\mathbf{e}}_{23}$	−${\mathbf{e}}_1$	${\mathbf{e}}_3$	−${\mathbf{e}}_{123}$	−${\mathbf{e}}_{13}$
${\mathbf{e}}_3$	${\mathbf{e}}_3$	−${\mathbf{e}}_{13}$	−${\mathbf{e}}_{23}$	1	${\mathbf{e}}_{123}$	−${\mathbf{e}}_2$	−${\mathbf{e}}_1$	${\mathbf{e}}_{12}$
${\mathbf{e}}_{12}$	${\mathbf{e}}_{12}$	−${\mathbf{e}}_2$	${\mathbf{e}}_1$	${\mathbf{e}}_{123}$	−1	${\mathbf{e}}_{23}$	−${\mathbf{e}}_{23}$	−${\mathbf{e}}_3$
${\mathbf{e}}_{23}$	${\mathbf{e}}_{23}$	−${\mathbf{e}}_{123}$	−${\mathbf{e}}_3$	${\mathbf{e}}_2$	−${\mathbf{e}}_{13}$	−1	${\mathbf{e}}_{12}$	−${\mathbf{e}}_1$
${\mathbf{e}}_{13}$	${\mathbf{e}}_{13}$	−${\mathbf{e}}_3$	−${\mathbf{e}}_{123}$	${\mathbf{e}}_1$	${\mathbf{e}}_{23}$	−${\mathbf{e}}_{12}$	−1	${\mathbf{e}}_2$
${\mathbf{e}}_{123}$	${\mathbf{e}}_{123}$	${\mathbf{e}}_{23}$	−${\mathbf{e}}_{13}$	${\mathbf{e}}_{12}$	−${\mathbf{e}}_3$	−${\mathbf{e}}_1$	${\mathbf{e}}_2$	−1

.

Referring to (2) and (3), a ${\mathbb{G}}_3$ matrix with m-row and n-column, noted ${\mathbb{G}}_3^{m\times n}$, is constructed as follows [20]

$$\begin{array}{cc}\hfill \mathbf{A}& ={\mathbf{A}}_0+{\mathbf{A}}_1{\mathbf{e}}_1+{\mathbf{A}}_2{\mathbf{e}}_2+{\mathbf{A}}_3{\mathbf{e}}_3+{\mathbf{A}}_4{\mathbf{e}}_{12}\hfill \\ & +{\mathbf{A}}_5{\mathbf{e}}_{23}+{\mathbf{A}}_6{\mathbf{e}}_{13}+{\mathbf{A}}_7{\mathbf{e}}_{123},\hfill \end{array}$$

(5)

where ${\mathbf{A}}_k$ for $k=1,2,3,\dots ,7$ are all $m\times n$ real number matrices. The transpose with reversion of $\mathbf{A}$ is denoted by ${\mathbf{A}}^H$

$$\begin{array}{cc}\hfill {\mathbf{A}}^H& ={\mathbf{A}}_0^T+{\mathbf{A}}_1^T{\mathbf{e}}_1+{\mathbf{A}}_2^T{\mathbf{e}}_2+{\mathbf{A}}_3^T{\mathbf{e}}_3-{\mathbf{A}}_4^T{\mathbf{e}}_{12}\hfill \\ & -{\mathbf{A}}_5^T{\mathbf{e}}_{23}-{\mathbf{A}}_6^T{\mathbf{e}}_{13}-{\mathbf{A}}_7^T{\mathbf{e}}_{123},\hfill \end{array}$$

(6)

where ${\mathbf{A}}_i^T$ for $k=1,2,3,\dots ,7$ denotes the transpose.

2.3. G-MODEL

A compact polarized GA model for the vector-sensor array was proposed in [21], named G-MODEL, which models the six-component outputs of a vector sensor holistically using a multi-vector in ${\mathbb{G}}_3$. Suppose there are K narrow-band, far-field and uncorrelated sources with wavelength $\lambda$ impinging on an array, which includes Q vector sensors. Define ${\theta }_k\in [0,2\pi )$, ${\varphi }_k\in [0,\pi )$, ${\gamma }_k\in [0,\pi /2)$ and ${\eta }_k\in [-\pi ,\pi )$ are the azimuth angle, elevation angle, polarization amplitude angle and phase difference angle of the $kth$ source, respectively.

Define ${\mathit{u}}_{\mathit{k}}=cos{\theta }_{k}sin{\varphi }_k{\mathit{e}}_{\mathbf{1}}+sin{\theta }_{k}sin{\varphi }_k{\mathit{e}}_{\mathbf{2}}+cos{\varphi }_k{\mathit{e}}_{\mathbf{3}}$ as the unit vector (see Figure 1) of the $kth$ source when it impinges on the sensor at the origin. ${\mathit{v}}_{k1}=-sin{\theta }_k{\mathit{e}}_{\mathbf{1}}+cos{\theta }_k{\mathit{e}}_{\mathbf{2}}$ and ${\mathit{v}}_{k2}=cos{\theta }_{k}cos{\varphi }_k{\mathit{e}}_{\mathbf{1}}+sin{\theta }_{k}cos{\varphi }_k{\mathit{e}}_{\mathbf{2}}-sin{\varphi }_k{\mathit{e}}_{\mathbf{3}}$ are unit multi-vectors. The position vector of the $qth$ sensor is ${\mathit{r}}_q=r_{q1}{\mathit{e}}_{\mathbf{1}}+r_{q2}{\mathit{e}}_{\mathbf{2}}+r_{q3}{\mathit{e}}_{\mathbf{3}}$. The output of the $qth$ vector sensor in the array is denoted by [21]

$$Y_{EH}^{\left(q\right)}\left(t\right)=\sum _{k=1}^K{X}_q\left({\theta }_k,{\varphi }_k\right)V_k{P}_k{S}_k\left(t\right)+N_{EH}^{\left(q\right)}\left(t\right),$$

(7)

where $X_q\left({\theta }_k,{\varphi }_k\right)=e^{{\mathit{e}}_{123}\frac{2\pi }{\lambda }(cos{\theta }_{k}sin{\varphi }_k{r}_{q1}+sin{\theta }_{k}sin{\varphi }_k{r}_{q2}+cos{\varphi }_k{r}_{q3})}$ is the spatial phase factor of the $kth$ source incident on the $qth$ vector sensor.

$$\begin{array}{cc}& V_k=\left(1+{\mathit{u}}_k\right)\left[{\mathit{v}}_{k1},-{\mathit{v}}_{k2}\right],\hfill \\ & P_k=\left[\begin{array}{c}cos{\gamma }_k\\ sin{\gamma }_k{e}^{{\mathit{e}}_{123}{\eta }_k}\end{array}\right],\hfill \\ & S_k\left(t\right)=\left|S_k\left(t\right)\right|exp\left[{\mathit{e}}_{123}\left(2\pi f_{k}t\right)\right].\hfill \end{array}$$

In next section, the GA-ESPRIT algorithm is deduced based on the G-MODEL.

3. Proposed Algorithm

The basic premise of the ESPRIT algorithm is that there are identical subarrays, the spacing between subarrays is known and the structure of subarrays is identical, which satisfies the rotational invariance in space [13]. Uniform linear arrays (ULAs) appear when it comes to one-dimensional DOA estimation using conventional ESPRIT [1,13]. Compared with ULAs, double parallel uniform linear arrays (DPULAs) can identify two-dimensional DOA because of the special construction, which consists of two parallel ULAs [36,37,38]; therefore, the algorithm discussed in this paper is based on DPULAs.

3.1. Complex Representation Matrix and Related Calculations

In view of the paucity of research on calculations with multi-vector, the Complex Representation Matrix (CRM) [20] is introduced because of the mature matrix theories. Consider a matrix $\mathbf{A}\in {\mathbb{G}}_3^{m\times n}$, the CRM is defined by $\Psi \left(\mathbf{A}\right)$

$$\Psi \left(\mathbf{A}\right)=\left(\begin{array}{cc}{\mathbf{A}}_0+{\mathbf{A}}_3+({\mathbf{A}}_7+{\mathbf{A}}_4){\mathbf{e}}_{123}& \hfill -{\mathbf{A}}_1+{\mathbf{A}}_6+({\mathbf{A}}_2-{\mathbf{A}}_5){\mathbf{e}}_{123}\\ -{\mathbf{A}}_1-{\mathbf{A}}_6-({\mathbf{A}}_2+{\mathbf{A}}_5){\mathbf{e}}_{123}& \hfill {\mathbf{A}}_0-{\mathbf{A}}_3+({\mathbf{A}}_7-{\mathbf{A}}_4){\mathbf{e}}_{123}\end{array}\right).$$

(8)

Let $\nu =\left(-{\mathbf{e}}_1+{\mathbf{e}}_{13}\right)/2\in {\mathbb{G}}_3$, and its reversion is $\tilde{\nu }=\left(-{\mathbf{e}}_1-{\mathbf{e}}_{13}\right)/2\in {\mathbb{G}}_3$. Then,

$${\nu }^2={\tilde{\nu }}^2=0\mathrm{and}\nu \tilde{\nu }+\tilde{\nu }\nu =1,$$

(9)

which imply $\nu \tilde{\nu }\nu =\nu$, $\tilde{\nu }\nu \tilde{\nu }=\tilde{\nu }$, ${\left(\nu \tilde{\nu }\right)}^2=\nu \tilde{\nu }$, ${\left(\tilde{\nu }\nu \right)}^2=\tilde{\nu }\nu$.

It immediately follows that, for every $\mathbf{A}\in {\mathbb{G}}_3$, we have

$$\mathbf{A}={\mathbf{E}}_{2m}\Psi \left(\mathbf{A}\right){\mathbf{E}}_{2n}^H,$$

(10)

$$\Psi \left(\mathbf{A}\right)={\mathbf{Q}}_{2m}\left[\begin{array}{cc}\mathbf{A}\hfill & \mathbf{0}\hfill \\ \mathbf{0}\hfill & \mathbf{A}\hfill \end{array}\right]{\mathbf{Q}}_{2n},$$

(11)

where in (10) and (11) we have

$${\mathbf{E}}_{2k}=\left[\nu \tilde{\nu }{\mathbf{I}}_k\tilde{\nu }{\mathbf{I}}_k\right]\in {\mathbb{G}}_3^{k\times 2k},$$

(12)

$${\mathbf{Q}}_{2k}=\left[\begin{array}{cc}\nu \tilde{\nu }{\mathbf{I}}_k& \tilde{\nu }{\mathbf{I}}_k\\ \nu {\mathbf{I}}_k& \tilde{\nu }\nu {\mathbf{I}}_k\end{array}\right]\in {\mathbb{G}}_3^{2k\times 2k}.$$

(13)

${\mathbf{I}}_k$ denotes the $k\times k$ identity matrix. It is not difficult to prove that

$${\mathbf{Q}}_{2k}={\mathbf{Q}}_{2k}^H={{\mathbf{Q}}_{2k}}^{-1},$$

(14a)

$$\Psi \left({\mathbf{A}}^H\right)={(\Psi (\mathbf{A}\left)\right)}^H,$$

(14b)

$$\Psi \left({\mathbf{A}}^{+}\right)={(\Psi \left(\mathbf{A}\right))}^{+},$$

(14c)

where $\{+\}$ denotes the pseudo-inverse. Referring to (10) and (14c), the pseudo-inverse of any $\mathbf{A}\in {\mathbb{G}}_3$ is

$${\mathbf{A}}^{+}={\mathbf{E}}_{2n}{(\Psi \left(\mathbf{A}\right))}^{+}{\mathbf{E}}_{2m}^H.$$

(15)

Since ${\mathbf{e}}_{123}^2=-1$ and ${\mathbf{e}}_{123}$ commutes with all elements in ${\mathbb{G}}_3$, one can identify it with the complex imaginary unit j [20], and so we can view $\Psi \left(\mathbf{A}\right)$ given in (8) as a complex matrix.

3.2. Model for DPULAs

Consider a DPULA with $2M+2$ sensors, as shown in Figure 2, in which d and M refer to the spacing between two adjacent sensors and the number of sensors in per subarray, respectively. The array is divided into three subarrays. The 1st to Mth sensors on the x-axis compose the first subarray, the 2nd to $(M+1)$th sensors form the second subarray and the $(M+2)$th to $(2M+1)$th that located on a straight line parallel to the x-axis make up the third subarray. The reason for the division can be found in Figure 3, that is, there are two unknown DOA parameters in the model, which need two rotational invariance relations.

Since the three subarrays have the same structure and the same number of sensor, each output of them has only one phase difference for the same signal. Signals received by subarray one, two and three are defined as ${\mathbf{Y}}_{EH}^1$, ${\mathbf{Y}}_{EH}^2$ and ${\mathbf{Y}}_{EH}^3$, respectively. According to the above array model, the outputs of the three subarrays at time t are as follows

$$\begin{array}{cc}\hfill {\mathbf{Y}}_{EH}^1\left(t\right)& =\mathbf{A}\mathbf{S}\left(t\right)+{\mathbf{N}}^1\left(t\right),\hfill \\ \hfill {\mathbf{Y}}_{EH}^2\left(t\right)& =\mathbf{A}\mathbf{F}\mathbf{S}\left(t\right)+{\mathbf{N}}^2\left(t\right),\hfill \\ \hfill {\mathbf{Y}}_{EH}^3\left(t\right)& =\mathbf{A}\mathbf{G}\mathbf{S}\left(t\right)+{\mathbf{N}}^3\left(t\right),\hfill \end{array}$$

(16)

where

$$\begin{array}{cc}\hfill & {\mathbf{Y}}_{EH}^1\left(t\right)={\left[Y_{EH}^{\left(1\right)}\left(t\right),\dots ,Y_{EH}^{\left(M\right)}\left(t\right)\right]}^T,\hfill \\ \hfill & {\mathbf{Y}}_{EH}^2\left(t\right)={\left[Y_{EH}^{\left(2\right)}\left(t\right),\dots ,Y_{EH}^{(M+1)}\left(t\right)\right]}^T,\hfill \\ \hfill & {\mathbf{Y}}_{EH}^3\left(t\right)={\left[Y_{EH}^{(M+2)}\left(t\right),\dots ,Y_{EH}^{(2M+1)}\left(t\right)\right]}^T,\hfill \\ \hfill \end{array}$$

(17)

and

$$\begin{array}{cc}\hfill & \mathbf{A}=\left[a{\Gamma }_1,\dots ,a{\Gamma }_K\right],\hfill \\ \hfill & a\left({\Gamma }_{\mathit{k}}\right)={\left[1,\mathrm{x}\left({\theta }_k,{\varphi }_k\right),\dots ,x^{M-1}\left({\theta }_k,{\varphi }_k\right)\right]}^T{V}_k{P}_k,\hfill \\ \hfill & \mathrm{x}\left({\theta }_k,{\varphi }_k\right)=e^{{\mathit{e}}_{123}\frac{2\pi }{\lambda }\mathrm{d}cos{\theta }_{k}sin{\varphi }_k},\hfill \\ \hfill & \mathbf{F}=diag\left(f_1,\dots ,f_K\right),\mathbf{G}=diag\left(g_1,\dots ,g_K\right).\hfill \end{array}$$

(18)

According to (18), we find that the DOA information is contained in matrix $\mathbf{A}$, $\mathbf{F}$ and $\mathbf{G}$. Because $\mathbf{F}$ and $\mathbf{G}$ are diagonal matrices that only contain direction information of incident signals, the focus is the two matrices, i.e.,

$$\begin{array}{c}{f}_k=e^{{\mathit{e}}_{123}\frac{2\pi }{\lambda }\mathrm{d}cos{\theta }_{k}sin{\varphi }_k},\hfill \\ g_k=e^{{\mathit{e}}_{123}\frac{2\pi }{\lambda }\mathrm{d}sin{\theta }_{k}sin{\varphi }_k}.\hfill \end{array}$$

(19)

Clearly, it is easy to figure out the DOA in the light of (19) if we obtain the two ideal matrices $\mathbf{F}$ and $\mathbf{G}$. From the rules of subarray division, we can see that the latter $(M-1)$ sensors of subarray one and the former $(M+1)$ sensors of subarray two are overlapped. Thus, in order to reduce the computational complexity, subarray one and subarray two can be merged to form a new matrix ${\mathbf{Y}}_{EH}$, that is,

$${\mathbf{Y}}_{EH}\left(t\right)={\left[y_1\left(t\right),y_2\left(t\right),\dots ,y_{M+1}\left(t\right)\right]}^T.$$

(20)

After merging, the $(2M+2)th$ redundant sensor is added to subarray three to form a new subarray ${\mathbf{P}}_{EH}$, so that the third subarray has the same dimension as ${\mathbf{Y}}_{EH}$

$${\mathbf{P}}_{EH}\left(t\right)={\left[y_{M+2}\left(t\right),y_{M+3}\left(t\right),\dots ,y_{2M+2}\left(t\right)\right]}^T.$$

(21)

Let $\overline{\mathbf{A}}$ be the array flow pattern of ${\mathbf{Y}}_{EH}$, then

$$\begin{array}{c}\overline{\mathbf{A}}=\left[\overline{a}\right({\Gamma }_1),\overline{a}({\Gamma }_2),\dots ,\overline{a}({\Gamma }_K\left)\right],\hfill \\ \overline{a}\left({\Gamma }_k\right)={\left[1,x\left({\theta }_k,{\varphi }_k\right),\dots ,x^M\left({\theta }_k,{\varphi }_k\right)\right]}^T{V}_k{P}_k.\hfill \end{array}$$

(22)

${\mathbf{Y}}_{EH}$ and ${\mathbf{P}}_{EH}$ can be written as

$$\begin{array}{cc}\hfill {\mathbf{Y}}_{EH}\left(t\right)& =\overline{\mathbf{A}}\mathbf{S}\left(t\right)+{\mathbf{N}}_a\left(t\right),\hfill \\ \hfill {\mathbf{P}}_{EH}\left(t\right)& =\overline{\mathbf{A}}\mathbf{G}\mathbf{S}\left(t\right)+{\mathbf{N}}_b\left(t\right),\hfill \end{array}$$

(23)

where

$$\begin{array}{c}\hfill {\mathbf{N}}_a\left(t\right)=\left[\begin{array}{c}{\mathbf{N}}^1\left(t\right)\\ n_{M+1}\left(t\right)\end{array}\right],{\mathbf{N}}_b\left(t\right)=\left[\begin{array}{c}{\mathbf{N}}^3\left(t\right)\\ n_{2M+2}\left(t\right)\end{array}\right].\end{array}$$

Then, $\mathbf{B}\left(t\right)$ is defined as

$$\mathbf{B}\left(t\right)=\left[\begin{array}{c}{\mathbf{Y}}_{EH}\left(t\right)\hfill \\ {\mathbf{P}}_{EH}\left(t\right)\hfill \end{array}\right]=\mathbf{C}\mathbf{S}\left(t\right)+\mathbf{N}\left(t\right),$$

(24)

where

$$\begin{array}{c}\hfill \mathbf{C}=\left[\begin{array}{c}\overline{\mathbf{A}}\\ \overline{\mathbf{A}}\mathbf{G}\end{array}\right],\mathbf{N}\left(t\right)=\left[\begin{array}{c}{\mathbf{N}}_a\left(t\right)\\ {\mathbf{N}}_b\left(t\right)\end{array}\right].\end{array}$$

Finally, the output of the whole array is denoted by

$$\mathbf{B}\left(t\right)=\mathbf{C}\mathbf{S}\left(t\right)+\mathbf{N}\left(t\right).$$

(25)

3.3. Algorithm Details

It is assumed that the sources received by the vector-sensor array are random signals which are independent and uncorrelated. In the same way, the measuring noise on six antennas of each sensor is white noise with the same power.

3.3.1. Subspace Separation

Under the above assumption, theoretically, the covariance matrix of the array output is given by

$$\mathbf{R}=\mathbb{E}\left\{\mathbf{B}{\mathbf{B}}^H\right\}=\mathbf{C}{\mathbf{R}}_s{\mathbf{C}}^H+6{\sigma }^2{\mathbf{I}}_{2M+2},$$

(26)

where $\mathbb{E}\{·\}$ stands for the mathematical expectation operator, ${\sigma }^2$ is the noise power on each vector antenna, ${\mathbf{R}}_S=\mathbb{E}\left\{\mathbf{S}\left(t\right){\mathbf{S}}^H\left(t\right)\right\}$.

Since the geometric product is non-commutativity, the Eigenvalue Decomposition (ED) is different from the conventional real methods but similar to the quaternion case. In other words, there are two possible eigenvalues, namely the left and the right eigenvalue for ${\mathbb{G}}_3$ matrix. In the proposed algorithm, the right eigenvalue is selected because the right ED of ${\mathbb{G}}_3$ matrix can be converted to the right ED of its CRM [20].

The ED of $\mathbf{R}$ is denoted by

$$\mathbf{R}={\mathbf{U}}_s{\mathsf{\Sigma }}_s{\mathbf{U}}_s^H+{\mathbf{U}}_n{\mathsf{\Sigma }}_n{\mathbf{U}}_n^H.$$

(27)

According to the principle of subspace separation, ${\mathbf{U}}_s$ is the signal subspace corresponding to K larger eigenvalues, and ${\mathsf{\Sigma }}_s$ is a diagonal matrix composed of K larger eigenvalues. In addition, ${\mathbf{U}}_n$ is orthogonal to ${\mathbf{U}}_s$ and it is the noise subspace corresponding to the remaining $4(M+1)-K$ small eigenvalues. Similarly, ${\mathsf{\Sigma }}_n$ is a diagonal matrix composed of the remaining small eigenvalues.

In the actual processing, the received signal is usually sampled. So, for a certain number of snapshots N, (26) and (27) can be rewritten as

$$\begin{array}{cc}\hfill & \widehat{\mathbf{R}}=\frac{1}{N}\sum _{i=1}^N\mathbf{B}\left(t_i\right){\mathbf{B}}^H\left(t_i\right),\hfill \\ \hfill & \widehat{\mathbf{R}}={\widehat{\mathbf{U}}}_s{\widehat{\mathsf{\Sigma }}}_s{\widehat{\mathbf{U}}}_s^H+{\widehat{\mathbf{U}}}_n{\widehat{\mathsf{\Sigma }}}_n{\widehat{\mathbf{U}}}_n^H.\hfill \end{array}$$

(28)

Because the space formed by the eigenvectors corresponding to the larger eigenvalues is the same as the space formed by the steering multi-vectors of the incident signals, that is, $span\left\{{\mathbf{U}}_s\right\}=span\left\{\mathbf{C}\right\}$, there exists a unique non-singular matrix $\mathbf{T}$, which satisfies

$${\mathbf{U}}_s=\mathbf{C}\mathbf{T}.$$

(29)

The rotational invariance relations exist among three subarrays, but ${\mathbf{U}}_s$ is the signal subspace of the whole array; therefore, after obtaining ${\mathbf{U}}_s$, the signal subspace of three subarrays must be separated. By the arrangement of sensor array, we find that the signal subspace of three subarrays can be calculated by

$$\begin{array}{cc}\hfill {\mathbf{U}}_{s1}& ={\mathbf{K}}_1{\mathbf{U}}_s=\mathbf{C}\mathbf{T},\hfill \\ \hfill {\mathbf{U}}_{s2}& ={\mathbf{K}}_2{\mathbf{U}}_s=\mathbf{C}\mathbf{F}\mathbf{T},\hfill \\ \hfill {\mathbf{U}}_{s3}& ={\mathbf{K}}_3{\mathbf{U}}_s=\mathbf{C}\mathbf{G}\mathbf{T},\hfill \end{array}$$

(30)

where ${\mathbf{U}}_{s1}$, ${\mathbf{U}}_{s2}$ and ${\mathbf{U}}_{s3}$ are signal subspaces of subarray one, subarray two and subarray three, respectively.

$$\begin{array}{cc}\hfill {\mathbf{K}}_1& ={\left[\begin{array}{cc}{\mathbf{I}}_M\hfill & {\mathbf{0}}_{M\times (M+2)}\hfill \end{array}\right]}_{M\times (2M+2)},\hfill \\ \hfill {\mathbf{K}}_2& ={\left[\begin{array}{ccc}{\mathbf{0}}_{M\times 1}\hfill & {\mathbf{I}}_M\hfill & {\mathbf{0}}_{M\times (M+1)}\hfill \end{array}\right]}_{M\times (2M+2)},\hfill \\ \hfill {\mathbf{K}}_3& ={\left[\begin{array}{ccc}{\mathbf{0}}_{M\times (M+1)}\hfill & {\mathbf{I}}_M\hfill & {\mathbf{0}}_{M\times 1}\hfill \end{array}\right]}_{M\times (2M+2)}.\hfill \end{array}$$

(31)

3.3.2. Rotation Invariance

From (30), the pivotal matrices $\mathbf{F}$ and $\mathbf{G}$ can be found. So, let

$${\mathbf{U}}_{s2}={\mathbf{U}}_{s1}{\mathsf{\Psi }}_x$$

(32)

in the same way,

$${\mathbf{U}}_{s3}={\mathbf{U}}_{s1}{\mathsf{\Psi }}_y$$

(33)

It is discovered that the eigenvalues of ${\mathsf{\Psi }}_x$ and ${\mathsf{\Psi }}_y$ are diagonal elements of $\mathbf{F}$ and $\mathbf{G}$, respectively.

Equations (32) and (33) are equations themselves, and are usually solved by the Least Squares (LS) method [29,36,38,39]; however, LS only takes the error on the left side of the equation into account, it ignores that the coefficient matrix also has an error; therefore, in order to reduce the error caused by solving the equation as much as possible, this paper considers a more accurate method—TLS [13]. Next, the solution of the equation is obtained by taking (32) as an example.

Combining the idea of TLS with the orthogonal property of subspace, we define a new matrix ${\mathbf{U}}_{s12}=\left[{\mathbf{U}}_{s1}{\mathbf{U}}_{s2}\right]$. In fact, the main aim is to seek a unitary matrix $\mathbf{D}\in {\mathbb{G}}_3^{M\times 2K}$, which is orthogonal to ${\mathbf{U}}_{s12}$. In other words, the space formed by $\mathbf{D}$ is orthogonal to the space formed by the column vectors of ${\mathbf{U}}_{s1}$ or ${\mathbf{U}}_{s2}$. So the $\mathbf{D}$ can be obtained from the ED of ${\mathbf{U}}_{s12}^H{\mathbf{U}}_{s12}$ [40]

$${\mathbf{U}}_{s12}^H{\mathbf{U}}_{s12}=\mathbf{E}\mathsf{\Lambda }{\mathbf{E}}^H,$$

(34)

where $\mathsf{\Lambda }$ is the diagonal matrix whose diagonal elements are composed by K multi-vectors that only have 0-grade-vector (can regard as non-zero real number) and $3K$ multi-vectors that equal to 0. $\mathbf{E}$ can be written as

$$\mathbf{E}=\left[\begin{array}{cc}{\mathbf{E}}_{11}\hfill & {\mathbf{E}}_{12}\hfill \\ {\mathbf{E}}_{21}\hfill & {\mathbf{E}}_{22}\hfill \end{array}\right].$$

(35)

Let ${\mathbf{E}}_N=\left[\begin{array}{c}{\mathbf{E}}_{12}\hfill \\ {\mathbf{E}}_{22}\hfill \end{array}\right]$, which is composed by eigenvectors whose eigenvalues are 0 and form the noise subspace. Since ${\mathbf{U}}_{s12}$ is signal subspace, we find that $\mathbf{D}={\mathbf{E}}_N$, i.e.,

$${\mathbf{U}}_{s12}\mathbf{D}=\left[{\mathbf{U}}_{s1}{\mathbf{U}}_{s2}\right]\left[\begin{array}{c}{\mathbf{E}}_{12}\hfill \\ {\mathbf{E}}_{22}\hfill \end{array}\right]=0.$$

(36)

Then,

$${\mathsf{\Psi }}_x=-{\mathbf{E}}_{12}{\mathbf{E}}_{22}^{+}.$$

(37)

The pseudo-inverse of ${\mathbb{G}}_3$ matrix ${\mathbf{E}}_{22}$ can be found in (15).

3.3.3. Angle Estimation

The azimuth and elevation angle of K signals are included in $\mathbf{F}$ and $\mathbf{G}$. In theory, the eigenvectors obtained by ED of these two matrices are both $\mathbf{T}$; however, in the actual calculation process, the two eigenvalue decomposition operations are carried out independently, which can not ensure that the arrangement of eigenvectors in them is reflected well; therefore, the diagonal elements of $\mathbf{F}$ and $\mathbf{G}$ should be matched.

Suppose that $\mathbf{T}\mathbf{1}$ and $\mathbf{T}\mathbf{2}$ are eigenvector matrices derived from GA-ED of ${\mathsf{\Psi }}_x$ and ${\mathsf{\Psi }}_y$, respectively. Then

$$\mathbf{O}=|\mathbf{T}{\mathbf{2}}^H\mathbf{T}\mathbf{1}|$$

(38)

where $\left\{\right|·\left|\right\}$ is the operator that gets magnitude of every multi-vector in a matrix. For the same signal, the eigenvectors in $\mathbf{T}\mathbf{1}$ and $\mathbf{T}\mathbf{2}$ corresponding to matched $f_k$ and $g_k$ are related; therefore, the order of diagonal elements in $\mathbf{F}$ and $\mathbf{G}$ can be adjusted by the coordinate of the largest element in each row (or column) of $\mathbf{O}$ to complete matching.

After observing (19), f and g are multi-vectors that only have scalar and 3-grade-vector, if we replace ${\mathit{e}}_{123}$ with the imaginary unit j of complex number, f and g can be regarded as complex numbers. Finally, we calculate ${\theta }_k$ and ${\varphi }_k$ with $f_k$ and $g_k$, that is,

$$\begin{array}{c}{\theta }_k={tan}^{-1}\left[\frac{angle\left(g_k\right)}{angle\left(f_k\right)}\right],\hfill \\ {\varphi }_k={sin}^{-1}\left\{\frac{\lambda }{2\pi }sqrt\left[angle{\left(g_k\right)}^2+angle{\left(f_k\right)}^2\right]\right\},\hfill \end{array}$$

(39)

where angle$(·)$ is the operator for getting phase angle. In conclusion, the steps of the GA-ESPRIT algorithm are:

• The original data received from three subarrays are integrated into the measurement model of the whole array according to (25);
• Calculate the covariance matrix $\widehat{\mathbf{R}}$, and then the ED in GA of $\widehat{\mathbf{R}}$ is performed and the signal subspace ${\mathbf{U}}_s$ can be obtained by the larger eigenvalues;
• According to (30), the signal subspace ${\mathbf{U}}_s$ of the whole array is divided into three subspaces ${\mathbf{U}}_{s1}$, ${\mathbf{U}}_{s2}$ and ${\mathbf{U}}_{s3}$;
• ${\mathsf{\Psi }}_x$ and ${\mathsf{\Psi }}_y$ can be obtained using TLS in GA, and details can be found in (34)–(37);
• The ED of ${\mathsf{\Psi }}_x$ and ${\mathsf{\Psi }}_y$ is performed to obtain matrices $\mathbf{F}$ and $\mathbf{G}$;
• The eigenvalues are matched in line with (38) and then taken them into Equation (39) to calculate K pairs direction angles.

Further, the corresponding relationship between the logic flow and steps of GA-ESPRIT is shown in Figure 4.

3.4. Complexity Analysis

As discussed in [21,22,27], the estimation of the data covariance matrix is an important factor to illustrate the complexity of ESPRIT algorithm and another one is ED, because they imply many repetitive operations and results, which mean heavy computational burden and memory requirements. Thus, we evaluate the time complexity of the two processes and space complexity in terms of real value memory requirements.

Suppose that an array composed of M vector sensors, and N snapshots are taken. LV-ESPRIT [13] and GA-ESPRIT consider six-component measurements of each vector sensor, whereas Q-ESPRIT [25] only records two-component measurements (electric field on x-axis and y-axis.); therefore, we compare the complexity between LV-ESPRIT and GA-ESPRIT. The output of each vector sensor for each signal consists of six complex numbers in LV-ESPRIT, while GA-ESPRIT only has one multi-vector with vector and bivector parts.

The geometric product of two multi-vectors received by the array output implies 36 real multiplications [21], which is nine times as many real multiplications as two complex numbers. As mentioned in Section 2, the ED of a ${\mathbb{G}}_3$ matrix is calculated by its CRM; therefore, the time complexity of the two algorithms is shown in Table 1. As for space complexity, the memory requirements of a real number is used to measure [21]. In the following two tables, CM is the Covariance Matrix, R represents real number.

The complexity comparison of these two algorithms can be found in

Table 2. Complexity of GA-ESPRIT and LV-ESPRIT.

Method	Time Complexity		Space Complexity
	CM	ED	CM (R)	Eigenvalue (R)	Eigenvectors (R)
LV-ESPRIT	$O\left\{N·{\left(6M\right)}^2\right\}$	$O\left\{{\left(6M\right)}^3\right\}$	72$M^2$	6M	36$M^2$
GA-ESPRIT	$O\left\{9N·M^2\right\}$	$O\left\{{\left(2M\right)}^3\right\}$	8$M^2$	2M	64$M^2$

, where CM and R represent covariance matrix and real number, respectively. Observing the time complexity in Table 2, it is not difficult to find that the computational burdens of CM and ED in GA-ESPRIT are a quarter and 1/27 of these in LV-ESPRIT, respectively. As for space complexity, GA-ESPRIT achieves such a significant reduction, more than 1.5 times compared to LV-ESPRIT, which means that the memory pressure is alleviated, especially for the large data size. The reason for the above comparison results is the natural advantage of GA matrix operations. In detail, because the six-dimensional measurement data in LV-MODEL (stored as 12 real numbers) are mapped into a multi-vector in the G-MODEL (stored as six real numbers), the amount of calculation will be reduced to varying degrees with different matrix operations, which will also bring fewer data storage requirements. The superior description and calculation ability of GA for multi-dimensional signals make GA-ESPRIT a very notable method for direction finding.

4. Simulation Results and Analysis

In this section, we simulate and analyze the proposed GA-ESPRIT based on DPULAs with $d=\lambda /2$, discuss its feasibility and performance compared with LV-ESPRIT [14] (in complex number field) and Q-ESPRIT [24] (in quaternion field). The estimation accuracy is evaluated by Root Mean Square Error (RMSE), which is calculated by the average of 200 Monte Carlo simulation experiments.

The RMSE of DOA estimation is defined as

$$RMSE=\frac{1}{K}\sum _{k=1}^K\sqrt{\frac{1}{200}\sum _{k=1}^{200}{\left[(\Delta {\theta }_k\right]}^2+{\left(\Delta {\varphi }_k\right)}^2]},$$

(40)

where K, $\Delta {\theta }_k$ and $\Delta {\varphi }_k$ denote the number of incident signals and errors between the result calculated by DOA algorithm and direction angle initially defined in the experiment, respectively.

In actual applications, the sensor model errors [27,41,42,43] cannot be ignored, which main include sensor-position error, gain error and phase error. The sensor-position error, as defined in [27], is the error between the actual position and the ideal position of each vector sensor. In the simulation experiment, the sensor-position error is modeled as additive noise with uniform distribution in a certain range, that is,

$${\overline{\mathbf{k}}}_m={\mathbf{k}}_m+d\sqrt{P_{pe}}{\left[{\epsilon }_{mx},{\epsilon }_{my},{\epsilon }_{mz}\right]}^T$$

(41)

where ${\overline{\mathbf{k}}}_m$ and ${\mathbf{k}}_m$ are the actual position and ideal position of the $mth$ sensor in vector-sensor array, respectively. ${\epsilon }_{mx}$, ${\epsilon }_{my}$ and ${\epsilon }_{mz}$ are uniformly distributed noise terms. $P_{pe}$ represents the perturbation power of sensor-position error and the larger $P_{pe}$ means the greater deviation of the sensor from its ideal position. Further, referring to [43], the array output with the gain and phase error is denoted by

$$\mathbf{B}\left(t\right)=(\mathbf{I}+\mathsf{\Pi }\mathsf{\Xi })\mathbf{C}\mathbf{S}\left(t\right)+\mathbf{N}\left(t\right),$$

(42)

where

$$\begin{array}{c}\mathsf{\Pi }=diag\left({\eta }_1,{\eta }_2,\dots ,{\eta }_{2M+2}\right),\hfill \\ \mathsf{\Xi }=diag\left(exp\left({\mathit{e}}_{123}{\xi }_1\right),\dots ,exp\left({\mathit{e}}_{123}{\xi }_{2M+2}\right)\right),\hfill \end{array}$$

in which ${\eta }_i$ and ${\xi }_i$ ($i=1,2,3,\dots ,2M+2$) are gain error and phase disturbance, respectively. In this paper, we also model them as additive noise. In addition, the six components of all EMVSs are added with noise according to the Signal-to-Noise ratio (SNR) in the following experiments. The SNR is defined as $SNR=10lg(P_s/P_n)$, in which $P_s$ and $P_n$ are the power of signal and noise on each component, respectively.

In the first experiment, we consider three far-field, narrow-band and uncorrelated signals with parameters Γ $=\left\{{160}^{\circ },{80}^{\circ },{35}^{\circ },-{60}^{\circ }\right\}$, $\left\{{60}^{\circ },{50}^{\circ },{35}^{\circ },{60}^{\circ }\right\}$ and $\{{20}^{\circ },{110}^{\circ }$, ${45}^{\circ },{80}^{\circ }\}$ with respect to Signal-to-Noise ratio (SNR) vary −10 dB to 20 dB in two different cases. In addition, we set $M=7$ and the snapshot number is 200. The aim of the first experiment was to examine the performance of GA-ESPRIT, LV-ESPRIT and Q-ESPRIT under different noise statistical characteristics. Figure 5a shows the estimation results of three algorithms when ideal Gaussian white noise is added, whereas, the noise in Figure 5b is related. It can be concluded that the three algorithms have very close accuracy of calculating DOA at high levels of SNR from Figure 5a,b, while with the lower SNR, GA-ESPRIT has higher accuracy over the other two and can achieve remove the correlation of noise partially.

In the second experiment, we compare the performance of GA-ESPRIT, LV-ESPRIT and Q-ESPRIT when the sensor-position error exists. Assume that two signals with Γ $=\left\{{58}^{\circ },{77}^{\circ },{35}^{\circ },-{60}^{\circ }\right\}$ and $\left\{{136}^{\circ },{50}^{\circ },{35}^{\circ },{60}^{\circ }\right\}$ impinge on a DPULA with $M=9$. Figure 6a shows the performance of the three algorithms when sensor-position error exists with different intensities. Meanwhile, we set SNR to 10 dB and the snapshot number is 200. The sensor-position error of the array sensor is changed by the value of $P_{pe}$, whose range is 0–0.07. It can be seen in Figure 6b that we fix $P_{pe}=0.02$ to observe the estimation of the three algorithms by altering SNR from −10 dB to 20 dB. Figure 6a,b both imply that accuracy of GA-ESPRIT is highest in the presence of the sensor-position error, so the conclusion is that GA-ESPRIT has the strongest robustness against sensor-position errors among the three algorithms.

The third experiment is also designed for two cases. Case one is that only gain error exists (see Figure 7a), while for case two, only phase error exists (see Figure 7b). Other conditions are the same as experiment two except that there is no position error. The gain error is constructed by the random numbers, whose mean value is 1 and variance is 0.2, and the phase error is constructed by the random numbers with zero-mean and 0.005 variance. We can learn from Figure 7a,b that, whether there is gain error or phase error, GA-ESPRIT can maintain the estimation accuracy very well, especially in low SNR.

In general, it is because Q-ESPRIT only takes part of the array output information into consideration that makes large RMSE. The reason for LV-ESPRIT’s poor accuracy in the face of the sensor-model error would be that its “long vector” destroys the orthogonality of the signal components. The improvement of detection robustness of GA-ESPRIT largely results from the fact that it can effectively preserve the orthogonality of the signal components and guarantee the completeness of the information.

5. Conclusions

In this paper, considering that the GA representation contains physical interpretations and complete information of incident signals, we use the idea of the traditional ESPRIT algorithm to find multiple EM signals in the direction finding method in GA. In particular, the model for DPULAs was built in GA and GA-ESPRIT was successfully derived using new calculation rules to achieve the two-dimensional DOA estimation. Compared with the previous ESPRIT algorithms, due to the robustness to sensor-model error and correlated noise, our proposed approach has potential in many practical situations, such as military radar in difficult environments. According to the experimental results, we have confirmed that the GA-ESPRIT has improved accuracy in two-dimension DOA estimation and can resist environmental interference to some extent. More importantly, the proposed algorithm achieves a reduction of more than 1/3 of the memory requirements while the time complexity is also greatly decreased.

Future works on the GA-ESPRIT will include polarization parameter estimation by optimizing matrix operations in GA and the ability of DOA recognition when facing coherent EM signals. It is expected that the proposed GA-ESPRIT will be an efficient DOA estimator.