Two-dimensional direction finding of coherent signals with a linear array of vector hydrophones

Abstract

In this paper, we propose a parallel factor (PARAFAC) analysis based two dimensional direction finding algorithm for coherent signals using a uniformly linear array of vector hydrophones. By forming a PARAFAC model using spatial signature of vector hydrophone array, the new algorithm requires neither spatial smoothing nor vector-field smoothing to decorrelate the signal coherency. We also establish that the azimuth-elevation angles of K coherent signals can be uniquely determined by PARAFAC analysis, as long as the number of hydrophones $L \ge 2K - 1$. In addition, because the vector hydrophone array manifold contains no phase information, this new algorithm can offer high azimuth-elevation estimation accuracy by setting vector hydrophones to space much farther apart than a half-wavelength. Simulation results are finally presented to verify the efficacy of the proposed algorithm.

1 Introduction

Estimation of azimuth-elevation angles using an array of vector hydrophones is of great importance in underwater acoustic signal processing and communications. A vector hydrophone consists of two or three orthogonally oriented velocity hydrophones plus a pressure hydrophone. Vector hydrophone measurement model was first introduced into the signal processing community in Nehorai and Paldi (1994). Since then, many efficient azimuth-elevation direction finding algorithms using vector hydrophone array have been presented (Wong and Zoltowski 1997a, b, 1999, 2000; Tichavsky et al. 2001; Zoltowski and Wong 2000; He and Liu 2009). These algorithms can achieve asymptotically efficient estimation performance in the case of incoherent signals. However, for real-world applications, the assumption of incoherent signals is often violated in scenarios where multipath exists. Coherent signals could reduce the rank of signal covariance matrix below the number of incident signals, and hence, degrade critically the algorithmic performance.

To deal with coherent signal direction finding, spatial smoothing technique (Shan et al. 1985; Wang and Liu 1998) is often used. In spatial smoothing, sensors are grouped into several (possibly overlapping) subarrays, all with the same subarray response. The data correlation matrices from these subarrays are averaged to decorrelate the coherency of the incoming signals. Preprocessing with spatial smoothing may reduce the effective array aperture, thereby degrading the array’s angular resolution and accuracy. For vector hydrophone array, a more sophisticated decoherency method is to perform smoothing across the particle-velocity-component’s data correlation matrices (Tao et al. 2008, 2007; He et al. 2010). However, these methods require a planar array geometry and/or 2D iteratively searching to estimate 2D directions of incoming signals. Another decoherency method is so-called subarray averaging (He and Liu 2009). With subarray averaging, subspace-based method can be used to estimate the vector hydrophone array manifolds. However, this method require additional computations to pair the estimated vector hydrophone array components.

In this paper, we propose a new two dimensional direction finding algorithm for coherent signals using a uniformly linear array of vector hydrophones. The proposed algorithm formulates a PARAFAC model by using the spatial signature of vector hydrophone array to extract vector hydrophone array manifolds by PARAFAC analysis, without requiring to perform spatial smoothing or vector-field smoothing to decorrelate the signal coherency. We also establish that the 2D directions of K coherent signals can be uniquely determined by PARAFAC analysis, as long as the number of hydrophones $L \ge 2K - 1$. In addition, because the vector hydrophone array manifold does not contain time-delay phase factor, the proposed algorithm may not suffer direction cyclical ambiguity when the spacing between hydrophones extends beyond a half-wavelength. The inherent structure of vector hydrophone also allows the array aperture extension to offer enhanced angle estimation precision.

Throughout the paper, scalar quantities are denoted by regular lowercase letters. Lowercase bold type faces are used for vectors and regular uppercase letters for matrices. Superscripts T, H, $*$, and $\dag$ represent the transpose, conjugate transpose, complex conjugate and pseudo inverse, respectively. $\otimes$ symbolizes the Kronecker-product operator. $\mathbf{0}_{m,n}$ and $\mathbf{I}_m$ denote the $m \times n$ zero matrix and $m \times m$ identity matrix, respectively.

2 Mathematical Data Model

We shall adopt the data model used in He and Liu (2009) for an array of L vector hydrophones receiving K coherent source signals. Each vector hydrophone has three components, i.e., two velocity hydrophones and a pressure hydrophone, co-located in space. The $3L \times 1$ data vector measured by the array at time t can be expressed as

$$\begin{aligned} \mathbf{z}(t)=\sum \limits _{k = 1}^{K} \mathbf{a}\left( \theta _k, \phi _k\right) s_k(t)+\mathbf{n} (t)=\mathbf{A} \mathbf{s}(t)+\mathbf{n}(t) \end{aligned}$$

(1)

where $\mathbf{n}(t)$ is $3L \times 1$ complex vector,

$$\begin{aligned} \mathbf{A}= & {} \left[ \mathbf{a}\left( \theta _1, \phi _1\right) , \cdots , \mathbf{a}\left( \theta _K, \phi _K\right) \right] \in ^{3L \times K}, \end{aligned}$$

(2)

$$\begin{aligned} \mathbf{a}\left( \theta _k, \phi _k\right)= & {} \mathbf{q}\left( \theta _k, \phi _k\right) \otimes \mathbf{c}\left( \theta _k, \phi _k\right) \in ^{3L \times 1} \end{aligned}$$

(3)

$$\begin{aligned} \mathbf{q}\left( \theta _k, \phi _k\right)= & {} \left[ 1, e^{-j\omega _k}, \cdots , e^{-j(L - 1)\omega _k}\right] ^T \in ^{L \times 1}, \end{aligned}$$

(4)

$$\begin{aligned} \omega _k= & {} \frac{2\pi }{\lambda }\left( d_x u_k + d_y v_k\right) , \end{aligned}$$

(5)

$$\begin{aligned} \mathbf{c}\left( \theta _k, \phi _k\right)= & {} \mathbf{c}_k = \left[ 1, u_k, v_k\right] ^T \in ^{3 \times 1}, \end{aligned}$$

(6)

$$\begin{aligned} \mathbf{s}(t)= & {} \left[ s_1(t), \cdots , s_K(t)\right] ^T = {\varvec{\beta }}s_1(t) \in ^{K \times 1} \end{aligned}$$

(7)

$$\begin{aligned} {\varvec{\beta }}= & {} \left[ \beta _1, \cdots , \beta _K\right] ^T \in ^{K \times 1} \end{aligned}$$

(8)

In above equations, $0 \le \theta _k < \pi /2$ denotes the $k\hbox {th}$ signal’s elevation angle, and $0 \le \phi _k < 2\pi$ represents the $k\hbox {th}$ signal’s azimuth angle. $u_k = \sin \theta _k \cos \phi _k$ and $v_k = \sin \theta _k \sin \phi _k$ signify the direction cosines along the x-axis and y-axis, respectively. $\mathbf{A}$ is the vector hydrophone array manifold, with $\mathbf{q}(\theta _k, \phi _k)$ being array spatial response vector, and $\mathbf{c}(\theta _k, \phi _k)$ being vector hydrophone’s response vector. $\beta _k$ is the complex attenuation coefficient of the $k\hbox {th}$ signal, with $\beta _k \ne 0$ and $\beta _1 = 1$. With a total of N data samples $\mathbf{z}(t_1), \cdots , \mathbf{z}(t_N)$ taken at the distinct time indices $\{t_n, n = 1, \ldots , N\}$, our objective is to determine $\{\theta _k, \phi _k, k = 1, \ldots , K\}$ from these samples.

3 Proposed Solution

From (1) and (7), $\mathbf{z}(t)$ can be reexpressed as

$$\begin{aligned} \mathbf{z}(t)= \mathbf{A} {\varvec{\beta }}s_1(t)+\mathbf{n}(t) = \mathbf{b} s_1(t) +\mathbf{n}(t) \end{aligned}$$

(9)

where

$$\begin{aligned} \mathbf{b} = \sum _{k = 1}^K \beta _k \mathbf{a}(\theta _k, \phi _k) \in ^{3L \times 1} \end{aligned}$$

(10)

is spatial signature of vector hydrophone array, which contains sufficient information on signal direction $(\theta _k, \phi _k)$. The proposed algorithm is based on the formulation of a PARAFAC model from $\mathbf{b}$.

3.1 PARAFAC Model Formulation

The PARAFAC analysis is a multi-way method originating from psychometrics (Harshman 1970; Carroll and Chang 1970). PARAFAC model, as a useful data analysis tool, is a generalization of low-rank matrix decomposition to three-way arrays (TWAs) or multi-way arrays (MWAs). During the past decade, the PARAFAC model is gaining more and more interest in numerous and diverse applications, such as in sensor array signal processing (Sidiropoulos et al. 2000) and communications (Sidiropoulos et al. 2000). To formulate the PARAFAC model, we need reshape the $3L \times 1$ vector $\mathbf{b}$ to a $3 \times L$ matrix as

$$\begin{aligned} \mathbf{Z} = \underbrace{\left[ \mathbf{c}\left( \theta _1, \phi _1\right) , \cdots , \mathbf{c}\left( \theta _K, \phi _K\right) \right] }_{\mathbf{C}} \underbrace{\left[ \beta _1 \mathbf{q}\left( \theta _1, \phi _1\right) , \cdots , \beta _K \mathbf{q}\left( \theta _K, \phi _K\right) \right] ^T}_{\mathbf{Q}} = \mathbf{CQ}^T \end{aligned}$$

(11)

From (11), with a total of L vector hydrophones, we can form $P (P \ge 2)$ different spatial-shifted data sets, where each associates with $(L - P + 1)$ vector hydrophones. The $p\hbox {th}$ spatial-shifted data set has the form

$$\begin{aligned} \mathbf{Z}_p = \mathbf{C}{\varvec{\Phi }}_p \mathbf{Q}_1^T, \quad p = 1, \ldots , P \end{aligned}$$

(12)

where

$$\begin{aligned} {\varvec{\Phi }}_p = \mathrm{diag}\left[ e^{-j \omega _1 (p - 1)}, \ldots , e^{-j \omega _K (p - 1)}\right] \end{aligned}$$

(13)

is a diagonal matrix and $\mathbf{Q}_1$ is the first $(L - P + 1)$ rows of $\mathbf{Q}$. Then, for $p = 1, \ldots , P$, we will have P different data sets $\{\mathbf{Z}_1, \ldots , \mathbf{Z}_P\}$. Note that these P data sets differ from each other only because the matrices ${\varvec{\Phi }}_p$ differ from one set to another. Next, let us stack these P matrices to form a three-way array $\underline{\mathbf{Z}}$, of which the $(i,j,p)\hbox {th}$ element can be written as

$$\begin{aligned} z_{i,j,p} = [\underline{\mathbf{Z}}]_{i,j,p} = \sum _{k = 1}^{K} c_{i,k} {\varphi }_{p,k} \alpha _{j,k} \end{aligned}$$

(14)

where $i = 1, \ldots , 3$, $j = 1, \ldots , N-P+1$, $p = 1, \ldots , P$, $c_{i,k}$ and $\alpha _{j,k}$, respectively, denote the $(i,k)\hbox {th}$ and the $(j,k)\hbox {th}$ entries of $\mathbf{C}$ and $\mathbf{Q}_1$, $\varphi _{p,k}$ represents the $(p,k)\hbox {th}$ element of the matrix ${\varvec{\Psi }}$ defined as

$$\begin{aligned} {\varvec{\Psi }} = \left[ \begin{array}{cccc} 1 &{}\quad 1 &{}\quad \cdots &{}\quad 1 \\ e^{-j \omega _1} &{}\quad e^{-j \omega _2} &{}\quad \ldots &{} e^{-j \omega _K} \\ \vdots &{}\quad \vdots &{}\quad \cdots &{}\quad \vdots \\ e^{-j \omega _1 (P-1)} &{}\quad e^{-j \omega _2 (P-1)} &{}\quad \ldots &{} e^{-j \omega _K (P-1)} \\ \end{array} \right] \end{aligned}$$

(15)

Apparently, the matrices ${\varvec{\Phi }}_p$ and ${\varvec{\Psi }}$ have the following relationship

$$\begin{aligned} {\varvec{\Phi }}_p = \mathcal{D}_p\{{\varvec{\Psi }}\}, \ \forall p = 1, \ldots , P \end{aligned}$$

(16)

where $\mathcal{D}_p\{ \cdot \}$ denotes the operator which takes the $p\hbox {th}$ row of the matrix in brackets and produces a diagonal matrix by placing this row on the main diagonal.

Equation (14) expresses $z_{i,j,p}$ as a sum of K rank-1 triple products. If certain conditions are satisfied, Eq. (14) represents a unique low-rank decomposition of $\underline{\mathbf{Z}}$. Therefore, the direction finding problem can be reformulated as the problem of low-rank decomposition of the three-way array $\underline{\mathbf{Z}}$, which can be solved by PARAFAC fitting.

3.2 PARAFAC Model Identifiability

In this subsection, we discuss the conditions under which the trilinear decomposition of $\underline{\mathbf{Z}}$ is unique. The identifiability conditions to be discussed are based on the notion of the Kruskal rank of a matrix (Kruskal 1977).

Definition

The Kurskal rank (or k-rank) of a matrix $\mathbf{A}$ is $k_\mathbf{A}$, if and only if every $k_\mathbf{A}$ columns of $\mathbf{A}$ are linearly independent and either $\mathbf{A}$ has $k_\mathbf{A}$ columns or $\mathbf{A}$ contains a set of $k_\mathbf{A}+1$ linearly dependent columns. Note that the Kruskal rank is always less than or equal to the conventional matrix rank. If $\mathbf{A}$ is of full column rank, then it is also of full k-rank.

To establish the identifiability, we rewrite the three-way array $\underline{\mathbf{Z}}$ compactly via Khatri-Rao matrix product format by using the relationship in (16)

$$\begin{aligned} \bar{\mathbf{Z}} = \left[ \begin{array}{c} \mathbf{Z}_1 \\ \mathbf{Z}_2 \\ \vdots \\ \mathbf{Z}_P \\ \end{array} \right] = \left[ \begin{array}{c} \mathbf{C} \mathcal{D}_1({\varvec{\Psi }}) \\ \mathbf{C} \mathcal{D}_2({\varvec{\Psi }}) \\ \vdots \\ \mathbf{C} \mathcal{D}_P({\varvec{\Psi }}) \\ \end{array} \right] \mathbf{Q}_1^T = ({\varvec{\Psi }} \odot \mathbf{C}) \mathbf{Q}_1^T \end{aligned}$$

(17)

The identifiability results are based on the the following theorem.

Theorem 1

For an $I \times J \times K$ TWA $\underline{\mathbf{X}}$ with a typical element $x_{i,j,k} = \sum _{f = 1}^F a_{i,f} b_{j,f} c_{k,f}$, $i = 1, \ldots , I, \ j = 1, \ldots , J, \ k = 1, \ldots , K$, where $a_{i,f}$, $b_{j,f}$ and $c_{k,f}$ respectively stand for the $(i, f)\hbox {th}$, $(j, f)\hbox {th}$ and $(k, f)\hbox {th}$ elements of the $I \times F$, $J \times F$ and $K \times F$ matrices $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$. If for $F > 1$,

$$\begin{aligned} k_\mathbf{A} + k_\mathbf{B} + k_\mathbf{C} \ge 2F+2 \end{aligned}$$

(18)

then $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$ are unique up to inherently unresolvable permutation and scaling of columns (Sidiropoulos et al. 2000).

On the basis of Theorem 1, we have the following result, which gives sufficient conditions for the number of the data samples to guarantee the models in (17) to be identifiable.

Theorem 2

The sufficient identifiable conditions for the TWAs (17) is, $L \ge 2K-1$.

Proof

From the analysis given in the “Appendix”, we know that for an vector hydrophone, every two vector hydrophone response vectors with distinct azimuth-elevation angles are linear independent. This means that the inequality $k_\mathbf{C} \ge 2$ holds. Next, for uniformly linear array, matrices ${\varvec{\Psi }}$ and $\mathbf{Q}_1$ are of Vandermonde structures, and hence, they are of full column-rank. Therefore, the following two equalities hold

$$\begin{aligned} k_{\varvec{\Psi }}= & {} \mathrm{rank}\{{\varvec{\Psi }}\} = \min (K, P) \end{aligned}$$

(19)

$$\begin{aligned} k_{\mathbf{Q}_1}= & {} \mathrm{rank}\{\mathbf{Q}_1\} = \min (K, L-P+1) \end{aligned}$$

(20)

Substituting (19) and (20) into (18), we have the identifiable condition is

$$\begin{aligned} k_{\mathbf{C}} + \min (K, P) + \min (K, L-P+1) \ge 2K+2 \end{aligned}$$

(21)

We discuss the following four cases

• 1. $P < K$, $L - P + 1 < K$. This implies that $L < 2K - 1$. In this case, in order to satisfy (21), we have $L \ge 2K-1$. This is contradictory with $L < 2K - 1$. Therefore, the TWA (17) is not identifiable for this case.

• 2. $P \ge K$, $L - P + 1 < K$. In this case, in order to satisfy (21), we have $L - P + 1 \ge K$. This is contradictory with $L -P + 1 < K$. Therefore, the TWA (17) is not identifiable for this case.

• 3. $P < K$, $L - P + 1 \ge K$. In this case, in order to satisfy (21), we have $P \ge K$. This is contradictory with $P < K$. Therefore, the TWA (17) is not identifiable for this case.

• 4. $P \ge K$, $L - P + 1 \ge K$. In this case, (21) becomes $k_\mathbf{C} + K + K \ge 2K+2$. Since $k_\mathbf{C} \ge 2$ is always satisfied, the TWA (17) is always identifiable for this case. Furthermore, $P \ge K$ and $L - P + 1 \ge K$ together lead to $L \ge 2K - 1$.

With the above discussions, the results in Theorem 2 are established. $\square$

3.3 PARAFAC Fitting

There exist several effective algorithms that can be used to fit the PARAFAC model. In this paper, the trilinear alternating least square (TALS) algorithm (Sidiropoulos et al. 2000) is adopted. In TALS, the parameters to be determined are separated into three sets, and each time a least squares cost function that depends only on one set is minimized. With the solution of this linear least squares problem, the subsequent stages of the TALS consist of applying the same principle on the remaining two sets of parameters. The TALS algorithm iterates, changing from one set to the next, until the variation of the loss function or of the parameters is less than a predefined convergence criterion. Since all the steps are optimizations in the least squares sense, the TALS algorithm is guaranteed to converge monotonically. In the simulations, the COMFAC algorithm (Sidiropoulos et al. 2000) is used to achieve a fast convergence of the TALS algorithm.

With the estimation of $\hat{\mathbf{C}}$ by PARAFAC fitting, the azimuth and elevation angle of the incoming signals can be easily obtained as follows. First, let the $k\hbox {th}$ column of $\hat{\mathbf{C}}$ be $\hat{\mathbf{c}}_k = [c_{k,1}, c_{k,2}, c_{k,3}]^T$. Then, the direction cosines $u_k$ and $v_k$ can, respectively, be estimated as $\hat{u}_k = c_{k,2}/c_{k,1}$ and $\hat{v}_k = c_{k,3}/c_{k,1}$. Finally, the elevation and azimuth angle estimates of the $k\hbox {th}$ signal are calculated as $\hat{\theta }_k = \sqrt{\hat{u}_k^2 + \hat{v}_k^2}$ and $\hat{\phi }_k = \angle (\hat{u}_k + \hat{v}_k)$.

It should be noted that the estimate $\hat{\mathbf{C}}$ is unique but up to some unknown permutation and scaling ambiguities. The scaling ambiguity can be easily resolved by normalizing each estimated column of $\hat{\mathbf{C}}$ with respect to its first element. The permutation ambiguity is unresolvable, but it is usually immaterial since the ordering of the estimated signal angles is unimportant.

3.4 Estimation of Spatial Signature of Vector Hydrophone Array

From the foregoing analysis, we can infer that with the estimation of spatial signature of vector hydrophone array $\hat{\mathbf{b}}$, the angles of signals can be estimated by forming a PARAFAC model from $\hat{\mathbf{b}}$, and then solving the formulated PARAFAC decomposition problem. Assuming temporally and spatially uncorrelated white complex Gaussian noise with zero-mean and variance $\sigma ^2$, the array covariance matrix can be expressed as

$$\begin{aligned} \mathbf{R} = E\{\mathbf{z}(t)\mathbf{z}^H(t)\} = \underbrace{E\{s_1(t)s_1^*(t)\}}_{\sigma ^2_s} \mathbf{b}{} \mathbf{b}^H + \sigma ^2 \mathbf{I}_{3L} \end{aligned}$$

(22)

Obviously, the rank of noiseless covariance matrix $\mathbf{R} - \sigma ^2 \mathbf{I}_{3L}$ is clearly equal to 1 due to the coherency of signals, and its eigenvalue decomposition is given by

$$\begin{aligned} \mathbf{R} - \sigma ^2 \mathbf{I}_{3L} = \mathbf{U} {\varvec{\Gamma }} \mathbf{U}^H = \lambda _1 \mathbf{u}_1 \mathbf{u}_1^H \end{aligned}$$

(23)

where $\mathbf{U} = [\mathbf{u}_1, \ldots , \mathbf{u}_{3L}]$, ${\varvec{\Gamma }} = \mathrm{diag}(\lambda _1, \ldots , \lambda _{3L})$, $\lambda _i$ and $\mathbf{u}_i$, are eigenvalues and eigenvectors, with $\lambda _1 > \lambda _2 = \lambda _3 = \cdots , = \lambda _{3L} = 0$.

From (22) and (23), we can infer that the vector $\mathbf{b}$ can be estimated by principal eigenvector of $\mathbf{R} - \sigma ^2 \mathbf{I}_{3L}$. In practical, the noise variance $\sigma ^2$ can be estimated by the average of $3L - 1$ smallest eigenvalues of $\mathbf{R}$.

3.5 Case for Partially Correlated or Uncorrelated Signals

The algorithm developed in previous subsections is based on the assumption that the incoming signals are fully correlated. In this subsection, we will show that this algorithm is also suitable for the scenario where the impinging signals are partially correlated or uncorrelated. For partially correlated or uncorrelated signals, the signal covariance matrix does not drop rank. To illustrate the adaptability of the proposed algorithm for this case, we let $\mathbf{E}_s = [\mathbf{e}_1, \ldots , \mathbf{e}_K]$ be the $3L \times K$ signal subspace matrix, whose columns are the $3L \times 1$ eigenvectors associated with the K largest eigenvalues of $\mathbf{R}$. Since $\mathbf{E}_s$ and $\mathbf{A}$ span the same space, there exists a $K \times K$ non-singular matrix $\mathbf{W}$ such that $\mathbf{E}_s = \mathbf{A} \mathbf{W}$. Thus, the principal eigenvector $\mathbf{e}_1$ of $\mathbf{R}$ equals $\mathbf{e}_1 = \mathbf{A} \mathbf{w}_1 = \sum _{k = 1}^K w_{1,k} \mathbf{a}(\theta _k, \phi _k)$, where $\mathbf{w}_1 = [w_{1,1}, \ldots , w_{1,K}]^T$ is the first column of $\mathbf{W}$. We can see that for scenario of the partially correlated or uncorrelated signals, the principal eigenvector $\mathbf{e}_1$ also gives an estimate of the spatial signature of vector hydrophone array, and therefore, the proposed algorithm is applicable to handle the partially correlated or uncorrelated signals as well. We will present a simulation to verify this point in Sect. 4.

3.6 Remarks

The presented algorithm obtains azimuth-elevation angle estimates by performing PARAFAC fitting to a TWA formulated from the spatial signature of vector hydrophone array. The superiorities of the proposed algorithm over other competitive algorithms are summarized as:

• 1. It requires neither spatial smoothing nor vector-field smoothing to decorrelate the signal coherency.

• 2. It uses a uniformly linear array to obtain two dimensional angle estimation.

• 3. The estimated azimuth and elevation angles are automatically paired without performing any pair matching computions.

4 Simulations

In this section, performance of the proposed algorithm is evaluated by comparing to the spatial smoothing algorithm in Wang and Liu (1998), the particle-velocity-field smoothing (PVFS) algorithm in Tao et al. (2008) and subarray averaging algorithm in He and Liu (2009). For the proposed algorithm and subarray averaging algorithm, 9-element uniformly linear array of vector hydrophone with inter-sensor spacing $d_x = d_y = d$ is used. For the spatial smoothing and PVFS algorithms, an L-shape array geometry, with 9 vector hydrophones uniformly placed vector hydrophones along x-axis for estimating $u_k$ and 9 vector hydrophones uniformly placed along y-axis for estimating $v_k$ is considered. The additive noise is assumed to be spatial white complex Gaussian, and the SNR is defined relative to each signal. The number of snapshots for each trial is $N = 256$. The result in each of the examples below is obtained from 500 independent Monte-Carlo trials.

Fig. 1

Direction cosine estimates for five signals. $\{u_1, u_2, u_3, u_4, u_5\}$ = $\{0.1, 0.2, 0.4, 0.5, 0.6\}$, $\{v_1, v_2, v_3, v_4, v_5\}$ = $\{0.15, 0.25, 0.35, 0.55, 0.65\}$

In the first example, we consider the scenario of five coherent signals impinging upon the above mentioned ULA. The signals are assumed to have the following direction cosines: $\{u_1, u_2, u_3, u_4, u_5\}$ = $\{0.1, 0.2, 0.4, 0.5, 0.6\}$, $\{v_1, v_2, v_3, v_4, v_5\}$ = $\{0.15, 0.25, 0.35, 0.55, 0.65\}$. The SNR is set as 20 dB. Figure 1 presents the scatter diagram of the direction cosine estimates for the five signals. This figure shows that the proposed algorithm successfully resolves the five coherent signals with 9 sensors, as stated in Theorem 2: $L \ge 2K - 1$.

Fig. 2

RMSEs of $u_1$ and $v_1$ estimates versus SNR. $u_1 = 0.1$, $v_1 = 0.2$, $u_2 = 0.2$, $v_2 = 0.3$. $N = 256$. 500 independent trials

In the second example, we examine the performance the proposed algorithm for the case of two signals with direction cosines $\{u_1, u_2\}$ = $\{0.1, 0.2\}$ and $\{v_1, v_2\}$ = $\{0.2, 0.3\}$. The root mean squared error (RMSE) of $u_1$ and $v_1$ estimates with different sensor spacing as a function of SNR are shown in Fig. 2. It is seen from the figure that the estimation errors of the proposed algorithm decrease as the sensor spacing increases. This phenomenon is in agreement with the intuitive reasoning that the larger of array aperture, the better the estimation performance. For conventional pressure-hydrophone array, increasing sensor spacing beyond a half-wavelength may lead to a set of cyclically ambiguous angle estimates. However, the proposed algorithm suffers no estimation ambiguity even the sensor spacing exceeding the Nyquist half-wavelength upper limit. This fact is because the proposed algorithm yields the angle estimates from the response vectors of the vector hydrophones.

In the third example, we plot the RMSEs of $u_1$ and $v_1$ estimates for different algorithms as a function of SNR in Fig. 3, where Cramér-Rao Bound (CRB) is also plotted for comparison. The simulation conditions are similar to those in the second example, except that the sensor spacing is $d = 4 \lambda$ for the proposed algorithm. The curves in figures unanimously demonstrate that the proposed algorithm can offer performance superior to those of the spatial smoothing algorithm, PVFS algorithm and subarray averaging algorithm, in terms of lower estimation errors.

Fig. 3

RMSEs of $u_1$ and $v_1$ estimates versus SNR. $u_1 = 0.1$, $v_1 = 0.2$, $u_2 = 0.2$, $v_2 = 0.3$. $N = 256$. 500 independent trials

In the last example, the RMSEs of $u_1$ and $v_1$ estimates for the proposed algorithm against the correlation coefficient $\rho \in [0, 1]$, which is defined as $\rho = \frac{E\{s_1(t) s_2^*(t) \}}{\sqrt{E[|s_1(t)|^2]E[|s_2(t)|^2]}}$, are shown in Fig. 4. We can see from the figure that the RMSEs of the direction cosine estimates remain almost unchanged with the varying of correlation coefficient $\rho$. This phenomenon is in agreement with the analysis given in Sect. 3.5 that the proposed algorithm is also applicable to handle the partially correlated or uncorrelated signals.

Fig. 4

RMSEs of $u_1$ and $v_1$ estimates versus correlation coefficient $\rho$. $u_1 = 0.1$, $v_1 = 0.2$, $u_2 = 0.2$, $v_2 = 0.3$. $N = 256$. 500 independent trials

5 Conclusions

We have presented a new two dimensional direction finding algorithm for coherent signals using a uniformly linear array of vector hydrophones. The PARAFAC model of spatial signature of vector hydrophone array have been formulated. The resulting estimator for estimating the azimuth and elevation angles can be obtained by PARAFAC fitting. The new algorithm requires neither spatial smoothing nor vector-field smoothing to decorrelate the signal coherency. The identifiability condition for PARAFAC fitting has been established as well. Finally, simulation results showing the superiority of the proposed algorithm have been presented.

Appendix

In this appendix, we will show that for an vector hydrophone, every two vector hydrophone response vectors with distinct azimuth-elevation angles are linear independent. For the problem under consideration, the values of $\theta _k$ and $\phi _k$ are confined as $\theta _k \in [0, \pi /2)$ and $\phi _k \in [0, 2\pi )$. Since $\mathbf{c}_k = \mathbf{c}_{\ell }$ if $\theta _k = \theta _{\ell } = 0$, regardless of the values of $\phi _k$ and $\phi _{\ell }$. Thus, the azimuth-elevation angles of the $k\hbox {th}$ and the $\ell \hbox {th}$ signals are distinct if

$$\begin{aligned} \left\{ \begin{array}{ll} (\theta _k, \phi _k) \ne (\theta _{\ell }, \phi _{\ell }) &{} \mathrm{if} \ \theta _k \ne 0 \\ \theta _k \ne \theta _{\ell } &{} \mathrm{if} \ \theta _k = 0 \\ \end{array} \right. \end{aligned}$$

Next, consider the following cases

• 1. $\theta _1 = 0$, $\phi _2 \in \{\pi /2, 3\pi /2\}$. In this case, the determinant of the matrix

  $$\begin{aligned} \mathrm{det}\left[ \begin{array}{cc} 1 &{}\quad 1 \\ \sin \theta _1 \sin \phi _1 &{}\quad \sin \theta _2 \sin \phi _2 \\ \end{array} \right] \ne 0 \end{aligned}$$

  Therefore, $\mathbf{c}_1$ and $\mathbf{c}_2$ are linear independent.

• 2. $\theta _1 = 0$, $\phi _2 \in \{0, \pi \}$. In this case, the determinant of the matrix

  $$\begin{aligned} \mathrm{det}\left[ \begin{array}{cc} 1 &{}\quad 1 \\ \sin \theta _1 \cos \phi _1 &{}\quad \sin \theta _2 \cos \phi _2 \\ \end{array} \right] \ne 0 \end{aligned}$$

  Therefore, $\mathbf{c}_1$ and $\mathbf{c}_2$ are linear independent.

• 3. $\theta _1 \ne 0$, $\theta _2 \ne 0$, $\phi _2 \ne m\pi /2, m = 0, 1, 2, 3$, $\phi _2 = \phi _1 + n \pi , n = 0, 1$. In this case, the determinant of the matrix

  $$\begin{aligned} \mathrm{det}\left[ \begin{array}{cc} 1 &{}\quad 1 \\ \sin \theta _1 \cos \phi _1 &{}\quad \sin \theta _2 \cos \phi _2 \\ \end{array} \right] \ne 0 \end{aligned}$$

  Therefore, $\mathbf{c}_1$ and $\mathbf{c}_2$ are linear independent.

• 4. $\theta _1 \ne 0$, $\theta _2 \ne 0$, $\phi _2 \ne m\pi /2, m = 0, 1, 2, 3$, $\phi _2 \ne \phi _1 + \pi$. In this case, the determinant of the matrix

  $$\begin{aligned} \mathrm{det}\left[ \begin{array}{cc} \sin \theta _1 \cos \phi _1 &{}\quad \sin \theta _2 \cos \phi _2 \\ \sin \theta _1 \sin \phi _1 &{}\quad \sin \theta _2 \sin \phi _2 \\ \end{array} \right] = \sin \theta _1 \sin \theta _2 (\cos \phi _1 \sin \phi _2 - \sin \phi _1 \cos \phi _2) \ne 0 \end{aligned}$$

  Therefore, $\mathbf{c}_1$ and $\mathbf{c}_2$ are linear independent.

Finally, we combine all the above results and conclude that for all distinct $(\theta _1, \phi _1)$ and $(\theta _2, \phi _2)$, the vectors $\mathbf{c}_1$ and $\mathbf{c}_2$ are linear independent.