Direction of Arrival Estimation for Non-Coherent Sub-Arrays via Joint Sparse and Low-Rank Signal Recovery

Estimating the directions of arrival (DOAs) of radio frequency sources using an array of antenna elements is an important problem that attracts much interest in different disciplines, such as signal processing and communications [1]. In many practical applications, e.g. when the sources and/or receivers are moving, it is required to estimate the DOAs from a single realization (snapshot) of the measurements. Several existing strategies can handle the single snapshot scenario. While the classical methods are based on minimization (often greedy) of the negative log-likelihood function [2, 3, 4] and more recent methods are based on training deep neural networks [5], another popular approach is based on formulating the problem as a sparse signal reconstruction task, and estimating the DOAs from the peaks of the magnitude of the recovered high-dimensional signal [6, 7, 8].

When the antenna elements are coherent,¹¹1Coherent elements are elements that are connected to the same local oscillator. Therefore, ignoring the noise, the phase differences between their observations are only due to the difference in the propagation delays between the location of the sources and the location of the elements. the DOA estimation performance improves with the increase of the array aperture [9, 10]. However, maintaining the coherence of all the elements of a large array is very demanding, especially at high carrier frequency, as it requires expensive hardware and a complicated calibration process that is prone to errors. Therefore, a possible alternative, with lower hardware complexity and cost, is to split a large array into several smaller sub-arrays, such that the coherence of the elements is kept only within each sub-array, while the large aperture of the entire array is exploited using signal processing techniques [11].

Most of the methods that can be applied to estimate the DOAs for non-coherent sub-arrays are either based on non-coherent processing of the sub-arrays (handling the sub-arrays as if they observe different signals) [12, 13, 14], or require a large number of snapshots with the same phase offsets between the sub-arrays (note, though, that unsynchronized local oscillators lead to different phase offsets in each snapshot) [15, 16, 17, 18, 19]. Therefore, to advance the performance of DOA estimation from a single snapshot of non-coherent sub-arrays, the work in [11] has proposed an approximate maximum likelihood approach. However, this approach assumes that the number of sources is known. Moreover, it is based on non-convex optimization with respect to an optimization variable whose dimension is the number of sources. Therefore, it requires a multidimensional search, or, when solved using iterative schemes, it depends on the quality of the initialization.

In this paper, we propose a novel method to estimate the DOAs from a single snapshot of non-coherent sub-arrays. It is based on convex optimization, and thus does not rely on good initializations and can easily handle scenarios with a large (possibly unknown) number of sources. Borrowing ideas from sparsity-based DOA estimation for coherent arrays [7] and from low-rank-based blind deconvolution [20], we formulate the problem as the reconstruction of a joint sparse and low-rank matrix and solve its convex relaxation. Then, the DOAs can be estimated from the solution of the convex problem. Yet, we show how an improvement is obtained if instead of estimating the DOAs, this solution is used to estimate the sub-arrays’ phase shifts, create a “phase-corrected” observation vector and use it to estimate the DOAs via another final (plain, coherent) sparsity-based DOA estimation. Numerical experiments show that the proposed approach outperforms other spectrum-based strategies, such as a variant of MUSIC [21] and a sparsity-based method that adapts self-calibration algorithms [22, 23] to the considered problem.

Consider $Q$ far-field sources, located at unknown angles $\bm{\theta}\triangleq[\theta_{1},\ldots,\theta_{Q}]^{T}$ and transmitting unknown narrowband signals that impinge on an $M$-element receiving array. Let the array be composed of $L$ sub-arrays, where the $\ell$th sub-array consists of $M_{\ell}$ elements ($\sum\limits_{\ell=1}^{L}M_{\ell}=M$). We assume that each sub-array is coherent, while different sub-arrays are non-coherent, i.e. there is an unknown phase shift between sub-arrays. Furthermore, we assume that only a single snapshot is available.

The waveform observed by each of the $L$ sub-arrays can be formulated as

$\displaystyle\bm{x}_{\ell}=\mathrm{e}^{-j\phi_{\ell}}\tilde{\bm{A}}_{\ell}(\bm{\theta})\tilde{\bm{s}}+\bm{n}_{\ell},\,\,\,\,\ell=1,\ldots,L,$

(1)

where $\phi_{\ell}$ is the unknown phase in the $\ell$th sub-array, $\tilde{\bm{A}}_{\ell}(\bm{\theta})\triangleq[\mathbf{a}_{\ell}(\theta_{1}),\ldots,\mathbf{a}_{\ell}(\theta_{Q})]\in\mathbb{C}^{M_{\ell}\times Q}$ is the $\ell$th sub-array manifold, where each column $\mathbf{a}_{\ell}(\theta)$ is the $\ell$th sub-array response to a signal which arrives at angle $\theta$. The vector $\tilde{\bm{s}}\in\mathbb{C}^{Q\times 1}$ represents the unknown signals, and $\bm{n}_{\ell}\in\mathbb{C}^{M_{\ell}\times 1}$ represents white, zero-mean, circular complex Gaussian noise. To simplify the formulation, we assume that the noise variance $\sigma^{2}$ is known, and equal at all receivers (the extension to non equal variances is straightforward). Regarding the steering vector $\mathbf{a}_{\ell}(\theta)$, assuming that the axis origin (i.e. the reference point) is defined in the middle of the entire array, and that the source angle is measured with respect to the boresight, then

$\displaystyle[\mathbf{a}_{\ell}(\theta)]_{i}=g_{\ell,i}(\theta)\mathrm{e}^{j\frac{2\pi}{\lambda}(x_{\ell,i}\mathrm{sin}\theta+y_{\ell,i}\mathrm{cos}\theta)},$

(2)

where $\lambda$ is the signal wavelength, $(x_{\ell,i},y_{\ell,i})$ denotes the antenna location, and $g_{\ell,i}$ denotes the antenna element pattern, which is assumed to be known. For example, for linear arrays of omnidirectional sensors we have that $g_{\ell,i}(\theta)=1$ and $y_{\ell,i}=0$.

To conclude this section, the problem at hand is to estimate the DOAs $\bm{\theta}$ from the sub-arrays’ observations $\{\bm{x}_{\ell}\}$ given in (1), where $\{\phi_{\ell}\}$ and $\tilde{\bm{s}}$ are nuisance parameters.

We propose to look at the problem from the sparse signal reconstruction perspective [7]. Namely, we define a DOA grid of length $N_{\theta}\gg Q$ (containing different hypotheses of a scalar DOA) and for each sub-array create the measurement matrix $\bm{A}_{\ell}\triangleq[\mathbf{a}_{\ell}(\theta_{1}),\ldots,\mathbf{a}_{\ell}(\theta_{N_{\theta}})]\in\mathbb{C}^{M_{\ell}\times N_{\theta}}$. The observation model in (1) can now be reformulated as

$\displaystyle\bm{x}_{\ell}=\alpha_{\ell}^{*}\bm{A}_{\ell}\bm{s}+\bm{n}_{\ell},\,\,\,\,\ell=1,\ldots,L,$

(3)

where $\alpha_{\ell}\triangleq\mathrm{e}^{j\phi_{\ell}}$ (the superscript ^∗ denotes the complex conjugate), and $\bm{s}\in\mathbb{C}^{N_{\theta}\times 1}$ obeys $\|\bm{s}\|_{0}=Q$, where $\|\cdot\|_{0}$ is the $\ell_{0}$ pseudo-norm that counts the number of non-zero entries of the vector. Note that the complicated non-linear model in (1) is replaced with the model in (3) that is bilinear with respect to the unknown parameters $\bm{\alpha}\triangleq[\alpha_{1},\ldots,\alpha_{L}]^{T}$ and $\bm{s}$. This suggests that, similarly to the “lifting” strategy used in [20] for the blind deconvolution problem (which is also bilinear), we can express (3) as a linear model with respect to a rank-1 matrix $\bm{Z}\triangleq\bm{s}\bm{\alpha}^{H}\in\mathbb{C}^{N_{\theta}\times L}$. Indeed,

$\displaystyle\bm{x}_{\ell}=\bm{A}_{\ell}\bm{Z}[:,\ell]+\bm{n}_{\ell},\,\,\,\,\ell=1,\ldots,L,$

(4)

where $\bm{Z}[:,\ell]$ denotes the $\ell$th column of $\bm{Z}$. Note that $\bm{Z}$ is not only rank-1, but also “row-sparse”, i.e., only $Q$ of its rows are non-zero. Note also that, in practice, the number of sources, $Q$, may not be known in advance. Therefore, a natural non-convex optimization problem that arises from (4) is

		$\displaystyle\underset{\bm{Z}}{\textrm{min}}\,\,\,\|\bm{Z}\|_{0,2}+\mu\cdot\mathrm{rank}(\bm{Z})$
		$\displaystyle\mathrm{s.t.}\,\,\,\,\sum\limits_{\ell=1}^{L}\left\|\bm{x}_{\ell}-\bm{A}_{\ell}\bm{Z}[:,\ell]\right\|_{2}^{2}\leq CM\sigma^{2},$		(5)

where $\mu$ and $C$ are positive parameters, $\|\cdot\|_{2}$ stands for the Euclidean norm, and $\|\cdot\|_{0,2}$ is a pseudo-norm that counts the number of non-zero rows of the matrix (or equivalently, counts the number of rows whose Euclidean norm is non-zero). Note that we currently ignore the fact that $\{\alpha_{\ell}\}$ have unit magnitude. Yet, this information will be used later on.

Hypothetically, given a solution to (3), one can extract from it an estimator of the signal $\bm{s}$, and estimate the DOAs as the angles associated with indices of the peaks of the magnitude of the latter estimator. However, (3) is a non-convex problem (due to its objective) that is not tractable to solve. Therefore, to obtain a convex problem that can be efficiently solved we perform convex relaxation on (3). Following the common practice [24, 25], we replace $\mathrm{rank}(\bm{Z})$ with the nuclear norm $\|\bm{Z}\|_{*}$ that sums the singular values of the matrix, and replace $\ell_{0}$-type pseudo-norms with their $\ell_{1}$-norm counterparts, i.e., we replace $\|\bm{Z}\|_{0,2}$ with $\|\bm{Z}\|_{1,2}\triangleq\sum\limits_{n=1}^{N_{\theta}}\sqrt{\sum\limits_{\ell=1}^{L}|Z[n,\ell]|^{2}}$. Therefore, we get the following convex optimization problem

		$\displaystyle\underset{\bm{Z}}{\textrm{min}}\,\,\,\|\bm{Z}\|_{1,2}+\mu\|\bm{Z}\|_{*}$
		$\displaystyle\mathrm{s.t.}\,\,\,\,\sum\limits_{\ell=1}^{L}\left\|\bm{x}_{\ell}-\bm{A}_{\ell}\bm{Z}[:,\ell]\right\|_{2}^{2}\leq CM\sigma^{2}.$		(6)

This convex program can be directly solved by the CVX Matlab package [26] (that applies the SDPT3 solver [27]). Since typically $L$ (the number of sub-arrays) is small, this procedure has tractable runtime (it took only a few seconds in our experiments with Intel-i7 laptop).

Let us denote by $\hat{\bm{Z}}_{\mathrm{cvx}}$ the solution to (3). Although the nuclear norm promotes low-rankness, it does not strictly impose rank-1. To strictly impose it, we apply the singular value decomposition (SVD) on $\hat{\bm{Z}}_{\mathrm{cvx}}$. Namely, we decompose $\hat{\bm{Z}}_{\mathrm{cvx}}=\bm{U}\bm{\Lambda}\bm{V}^{H}$, where $\bm{U}\in\mathbb{C}^{N_{\theta}\times N_{\theta}}$ and $\bm{V}\in\mathbb{C}^{L\times L}$ are unitary matrices and $\bm{\Lambda}$ is a rectangular diagonal matrix. Now, the rank-1 approximation of $\hat{\bm{Z}}_{\mathrm{cvx}}$ is given by $\hat{\bm{Z}}=\Lambda[1,1]{\bm{U}[:,1]}\bm{V}[:,1]^{H}$, where $\hat{\bm{s}}=\Lambda[1,1]\bm{U}[:,1]$ and $\hat{\bm{\alpha}}=\bm{V}[:,1]$ can be used as proxies to $\bm{s}$ and $\bm{\alpha}$.

We turn to develop two DOA estimation strategies. The first strategy, which has already been mentioned above, is to estimate the DOAs as the peaks of the magnitude of the $N_{\theta}\times 1$ vector $\hat{\bm{s}}$. Yet, recall that until now we ignored the fact that the “true” $\{\alpha_{\ell}\}$ have unit magnitude (since $\alpha_{\ell}=\mathrm{e}^{j\phi_{\ell}}$). To exploit this information, we propose another DOA estimation strategy (which improves the results in our experiments).

First, we estimate (only) the phase shifts $\{\phi_{\ell}\}$ from $\hat{\bm{Z}}$ by

$\displaystyle\hat{\phi}_{\ell}=\angle\hat{\alpha}[\ell]=\angle V[\ell,1],\,\,\,\,\ell=1,\ldots,L,$

(7)

where $\angle x$ stands for the phase of the complex number $x$. Looking back on the model in (1), we use the estimated phase shifts to obtain the “phase-corrected” observations of the entire array, denoted by $\hat{\overline{\bm{x}}}$, that can be (approximately) modeled as

$\displaystyle\hat{\overline{\bm{x}}}\triangleq\begin{bmatrix}\mathrm{e}^{j\hat{\phi}_{1}}\bm{x}_{1}\\ \vdots\\ \mathrm{e}^{j\hat{\phi}_{L}}\bm{x}_{L}\end{bmatrix}\approx\tilde{\bm{A}}(\bm{\theta})\tilde{\bm{s}}+\overline{\bm{n}},$

(8)

where $\tilde{\bm{A}}(\bm{\theta})\triangleq\left[\tilde{\bm{A}}_{1}^{T}(\bm{\theta})\ldots\tilde{\bm{A}}_{L}^{T}(\bm{\theta})\right]^{T}\in\mathbb{C}^{M\times Q}$ and $\overline{\bm{n}}\in\mathbb{C}^{M}$ is a noise vector, whose statistics are not expected to dramatically change (compared to $\{\bm{n}_{\ell}\}$) due to the phase shifts. As a minor remark, note that only the differences between the phase shifts estimates $\{\hat{\phi}_{\ell}\}$ matter, as the mean of the estimated phases can be absorbed into the unknown vector $\tilde{\bm{s}}$.

At this point, any DOA estimation method for coherent array that can handle the single snapshot case can be applied on the “phase-corrected” $\hat{\overline{\bm{x}}}$, modeled in (8). We choose to use the popular sparsity-based technique from [7]. Using the same notations of $\bm{s}$ and $\bm{A}_{\ell}$ that are used in (3), we reformulate (8) as

$\displaystyle\hat{\overline{\bm{x}}}\approx\bm{A}\bm{s}+\overline{\bm{n}},$

(9)

where $\bm{A}\triangleq\left[\bm{A}_{1}^{T}\ldots\bm{A}_{L}^{T}\right]^{T}\in\mathbb{C}^{M\times N_{\theta}}$, and estimate $\bm{s}$ as the minimizer of the following convex program

		$\displaystyle\underset{\bm{s}}{\textrm{min}}\,\,\,\|\bm{s}\|_{1}$
		$\displaystyle\mathrm{s.t.}\,\,\,\,\left\|\hat{\overline{\bm{x}}}-\bm{A}\bm{s}\right\|_{2}^{2}\leq CM\sigma^{2},$		(10)

where $C$ is a positive parameter and $\|\cdot\|_{1}$ is the $\ell_{1}$-norm. Again, this problem can be solved fastly using CVX [26], and we denote by $\hat{\bm{s}}_{\mathrm{cvx}}$ this solution. Finally, as done in [7], we estimate the DOAs as the peaks of the magnitude of the $N_{\theta}\times 1$ vector $\hat{\bm{s}}_{\mathrm{cvx}}$.

We perform computer simulations to evaluate the performance of the two proposed DOA estimation methods. The strategy where we estimate the DOAs from $\hat{\bm{Z}}$ will be referred to as Proposed1. The strategy where we estimate the phase shifts from $\hat{\bm{Z}}$ and then apply (quasi-)coherent sparsity-based method will be referred to as Proposed2. In all the experiments we use the parameters $\mu=1$ and $C=2$ in (3).

We compare our strategies with two reference methods. The first reference method is a sparsity-based approach that adapts the self-calibration algorithms [22, 23] to the considered non-coherent problem. Specifically, this technique obtains a matrix $\tilde{\bm{Z}}\in\mathbb{C}^{N_{\theta}\times L}$ using a convex program that is similar to (3) except that its objective includes only $\|\bm{Z}\|_{1,2}$ (i.e., $\mu=0$ and $C=2$). Then, following [23] (that has presented better results than [22]), the rank-1 approximation of $\tilde{\bm{Z}}$ is obtained using SVD. The DOAs are estimated as the peaks of the magnitude of the $N_{\theta}\times 1$ most dominant left singular vector of $\tilde{\bm{Z}}$. This method will be referred to as SparsityOnly. Observe that the difference between our Proposed1 and SparsityOnly is the fact that we include the nuclear norm in the convex optimization.

The second reference method is based on MUSIC [21]. It treats the realizations of the $L$ sub-arrays as $L$ snapshots of a single sub-array. We use this method, which will be referred to as MUSIC, with spatial smoothing [28] and forward-backward smoothing [29] (within each sub-array) that empirically improved its performance. Applying this non-coherent MUSIC is possible since in our experiments the sub-arrays are uniform linear arrays (ULAs) with similar sensor configuration. Note that all the examined methods are based on estimating the DOAs as the peaks of $N_{\theta}\times 1$ “spectrums”, and that only MUSIC requires the number of sources, $Q$, for generating its spectrum.

We consider a ULA of $M=24$ elements, with half wavelength spacing, partitioned into $L=4$ sub-arrays of size $M_{\ell}=6$. We start with $Q=2$ sources located at angles $\{0^{\circ},15^{\circ}\}$ (the boresight angle of the entire array is $0^{\circ}$) that transmit complex Gaussian signals with equal power. The unknown phase of each sub-array is selected at random from a uniform distribution over $[0,2\pi]$. To obtain statistical results, we perform $N_{\text{exp}}=250$ Monte Carlo experiments at each SNR between 0 dB to 30 dB. In each experiment we draw the Gaussian noise $\{\bm{n}_{\ell}\}$, the signal $\tilde{\bm{s}}$, and the unknown phases $\{\phi_{\ell}\}$. For each method, the performance is measured by the root mean square error (RMSE) of the $Q$ highest peaks.

The RMSE results are presented in Fig. 2. It can be seen that Proposed1, which promotes low-rankness within the convex optimization, it better than SparsityOnly. The best accuracy is observed for our Proposed2. The better performance of Proposed2 over Proposed1 (especially at high SNR) demonstrates the advantage of using the solution of the joint sparse and low-rank optimization problem for estimating the phase shifts rather than the DOAs.

Spectrum-based methods are especially preferable over maximum likelihood techniques when the number of sources is large. Thus, we turn to examine a scenario with $Q=4$ sources, located at angles $\{-15^{\circ},0^{\circ},15^{\circ},30^{\circ}\}$. The rest of the configuration remains as before. The RMSE results are presented in Fig. 2, and the estimators’ spectrums for one realization at SNR of 20 dB are displayed in Fig. 3. Comparing the spectrums of Proposed1 and SparsityOnly in Fig. 3, we see that the nuclear norm that is used in Proposed1 reduces the sparsity but significantly improves the accuracy of the DOAs estimates (i.e., the location of the peaks). We further see from Fig. 3 that MUSIC’s peaks are quite accurate for this realization, though, one source is almost not identified. Such detection failures explain the bad results of MUSIC at low to medium SNR in Fig. 2. Clearly, the best accuracy is observed for our Proposed2.

Note that in Fig. 2 the advantage of the two proposed methods over the reference methods is larger than in Fig. 2. This can be explained by the property that having more sources improves the phase offsets estimation (see the analysis in [11]). Specifically, this property directly explains the improved results compared to MUSIC (which is a non-coherent processing method). As for the improved advantage over SparsityOnly, note that the phase shifts estimation relates to the nuclear norm term in (3). While this term is used for both Proposed1 and Proposed2, it is eliminated ($\mu=0$) for SparsityOnly.

We addressed the problem of estimating the DOAs of multiple sources from a single snapshot of non-coherent sub-arrays. We formulated this problem as the reconstruction of a joint sparse and low-rank matrix and solved its convex relaxation. Unlike methods that are based on maximum likelihood [11], our approach can easily handle scenarios with a large (and possibly unknown) number of sources. We showed two ways to estimate the DOAs. The first strategy, estimates the DOAs directly from the solution of the convex problem and outperforms an existing sparsity-based method that does not promote low-rankness within the optimization. The second strategy, which has shown even improved results, uses the solution of the convex problem for estimating the phase shifts rather than the DOAs. The estimated phase shifts are then used to “correct” the observations, which allows DOA estimation with plain methods of coherent arrays. We note that there is work that questions the benefit of combining nuclear and $\ell_{1}$ norms for exact reconstruction of an entire high-dimensional matrix [30]. However, our paper shows that when the goal is to estimate from the observations only a rather small set of fundamental parameters (e.g., DOAs and phase shifts) then this combination of norms can be beneficial.