Structural-Missing Tensor Completion for Robust DOA Estimation with Sensor Failure ^†

Abstract

Array sensor failure poses a serious challenge to robust direction-of-arrival (DOA) estimation in complicated environments. Although existing matrix completion methods can successfully recover the damaged signals of an impaired sensor array, they cannot preserve the multi-way signal characteristics as the dimension of arrays expands. In this paper, we propose a structural-missing tensor completion algorithm for robust DOA estimation with uniform rectangular array (URA), which exhibits a high robustness to non-ideal sensor failure conditions. Specifically, the signals received at the impaired URA are represented as a three-dimensional incomplete tensor, which contains whole fibers or slices of missing elements. Due to this structural-missing pattern, the conventional low-rank tensor completion becomes ineffective. To resolve this issue, a spatio-temporal dimension augmentation method is developed to transform the structural-missing tensor signal into a six-dimensional Hankel tensor with dispersed missing elements. The augmented Hankel tensor can then be completed with a low-rank regularization by solving a Hankel tensor nuclear norm minimization problem. As such, the inverse Hankelization on the completed Hankel tensor recovers the tensor signal of an unimpaired URA. Accordingly, a completed covariance tensor can be derived and decomposed for robust DOA estimation. Simulation results verify the effectiveness of the proposed algorithm.

1. Introduction

Direction-of-arrival (DOA) estimation using sensor arrays plays an important role in various applications such as radar, sonar, and industrial Internet of Things (IIoT) [1,2,3,4,5]. Existing methods include subspace-based methods [6,7], sparsity-based methods [8,9,10], beamforming-based methods [11,12], statistical reconstruction-based methods [13,14], and neural-network-based methods [15,16]. Furthermore, advanced DOA estimation methods have mitigated the multi-path effects in practical environments, based on either iterative implementation approach [17,18] or low-rank decomposition and sparse representation approach [19]. However, all the aforementioned methods typically work under the presumption that the deployed sensor array remains unimpaired. Nevertheless, in practical environments with complex electromagnetic propagation and harsh working conditions, the gain-phase deviation or position deviation of sensors commonly occurs [20,21,22]. In more severe cases, some sensors will be totally damaged and become nonfunctional, i.e., the sensor failure problem. The above-mentioned methods are not robust to the sensor failure scenario, and thus, the DOA estimation performance will dramatically drop.

To achieve robust DOA estimation with sensor failure, matrix completion methods are developed to recover a complete signal matrix for the impaired array [23,24,25]. In particular, the principles of recursive least squares [23], nuclear-norm-based convex relaxation [26], and reweighted ${\ell }_{2,1}$-norm [24] can be incorporated with a low-rank regularization for signal matrix completion [27], leading to an effective DOA estimation. Meanwhile, co-array processing is adopted to help reconstruct an unimpaired array with improved performance [28]. To demonstrate the lower bound on the estimation accuracy under sensor failure, the Cramér–Rao bound (CRB) in various sensor failure scenarios is analysed for regular tetrahedral arrays [29]. Furthermore, the weighted sum of CRB for independent sensor failure event has been proposed for scenario of randomly impaired sensors [30]. Nevertheless, as the dimensions of sensor arrays expand, the matrix-based signal modelling and processing reflect at most two-dimensional (2-D) information, while complicated multi-way signal characteristics are ignored when it comes to multi-dimensional array processing. This phenomenon results in performance deterioration especially under sensor failure.

Different from matrix, tensor is utilized to preserve the original high-level information of data with three-dimensional (3D) structure and beyond 3D. Meanwhile, numerous tensor decomposition techniques such as the canonical polyadic decomposition (CPD) [31] and Tucker decomposition [32] have been proposed to retrieve parameters from the tensor data [33,34]. Based on tensor models and tensor decompositions, many tensor-based DOA estimation methods are proposed to enhance the estimation performance for multi-dimensional arrays [35,36]. Specifically, the CPD is applied to multiple-invariance array signals for effective DOA estimation [37]. The tensor subspace method is developed by formulating high-order singular value decomposition (HOSVD)-based spectral function [38]. The coupled CPD is proposed to fully exploit the spatial relevance of segmented signals to enhance the DOA estimation accuracy [39]. Moreover, to improve computational efficiency while maintaining estimation performance, the tensor train decomposition is developed to replace the CPD [40]. To increase the degrees-of-freedom for DOA estimation, the strategies of cross-correlation tensor processing [41], co-array tensor processing [42], and co-array tensor completion [43] have been investigated. Targeting the problem of coherent source DOA estimation in a multi-path environment, the tensor spatial smoothing [44] and tensor reconstruction [45] approaches are further designed to obtain a high robustness to the coherency of signal statistics. Unfortunately, none of these methods have considered the problem of sensor failure for multi-dimensional arrays. Therefore, it is of great significance to develop an effective tensor processing strategy tailored for robust DOA estimation with sensor failure.

From the perspective of recovering damaged signals as the matrix completion methods do, a promising solution is to complete tensor signals corresponding to the impaired multi-dimensional array. With regard to the tensor signal model of an impaired array, the missing elements will concentrate on entire fibers or even slices. This structural-missing pattern (i.e., general cases of signal-missing pattern including the fiber-missing and slice-missing pattern [46]) prevents the conventional low-rank tensor completion to be directly implemented, because the effective low-rank regularization requires a randomly missing pattern to guarantee sufficient spatial relevance [47]. With respect to the need for completing structurally missing elements, the slice completion methods are proposed exploiting multi-way delay-embedding transformation and Tucker rank increment [48,49]. Nevertheless, the complexity of rank-increment is high, and the Tucker model is unsuitable for the array tensor signal, which is normally modelled in a canonical polyadic (CP) format due to the deterministic source number. Furthermore, a CP approximation method has been proposed to reconstruct functional magnetic resonance imaging (fMRI) data with regular under-sampling patterns including both fiber-missing and slice-missing [50,51]. However, such reconstruction changes the undamaged pixels in the fMRI tensor, whereas the observed elements in the tensor signal are expected to remain unchanged to avoid extra model deviations including unwanted noise and phase error. As such, all these tensor completion methods cannot be applied to the problem of robust DOA estimation with sensor failure. Finding a successful approach to completing the structural-missing tensor signal remains a challenging task.

In this paper, we propose a structural-missing tensor completion algorithm for robust DOA estimation with sensor failure. In particular, the signal received at a uniform rectangular array (URA) is represented as a tensor, which contains whole fibers or slices of missing elements when several sensors are impaired. To enforce effective low-rank regularization on the structural-missing tensor signal, a spatio-temporal dimension augmentation method is developed to transform the original tensor into a dimensional-augmented Hankel tensor, whose missing elements are sufficiently dispersed and dimensional information is augmented. On this basis, the dimensional-augmented Hankel tensor can be optimized through low-rank completion, which is characterized by a convex tensor nuclear norm minimization problem. After solving this dimensional-augmented Hankel tensor completion problem, a completed tensor signal corresponding to an unimpaired array is obtained from an inverse Hankelization operation. As such, the completed covariance tensor can be derived and decomposed via CPD to achieve closed-form DOA estimation. According to our simulation results, the proposed algorithm shows high robustness to different simulated sensor failure scenarios compared to conventional methods.

Although the proposed algorithm is designed for the impaired URA with a uniform inter-sensor spacing, it can be extended to deployments of sparse arrays such as coprime array [52]. On one hand, the physical sparse array can be considered as a special case of partially damaged URA. Hence, the proposed algorithm can be applied to the incomplete sparse array signals for a direct completion. Based on this, the subsequent signal statistics can be derived and used for DOA estimation. On the other hand, for the augmented discontinuous co-array of the coprime array [53,54], the proposed algorithm can also be extended to complete the co-array signal. In this case, by incorporating the proposed algorithm in the co-array domain, the virtual array aperture can be further expanded, contributing to enhanced DOA estimation performance.

The preliminary results of this work were presented in the conference paper [55]. In this journal paper, we present more details of the tensor signal model and the spatio-temporal dimension augmentation method. Moreover, to ensure more effective tensor signal recovery, the dimensional-augmented Hankel tensor is completed with a low-rank regularization in this paper instead of Tucker approximation as in the conference paper [55]. We demonstrate the advantage of using the low-rank regularization both analytically and numerically. We elaborate the procedures of Hankel tensor completion and completed covariance tensor decomposition. Furthermore, to validate the effectiveness of the proposed algorithm, we present new simulation results, including comparison of estimation accuracy, comparison of angular resolution, performance of source identifiability, comparison with other completion methods, and comparison of computational time.

The rest of the paper is organized as follows. In Section 2, we formulate the incomplete tensor signal model under sensor failure. In Section 3, we design the structural-missing tensor completion method, which includes the procedures of spatio-temporal dimension augmentation and Hankel tensor completion. In Section 4, we propose the robust DOA estimation method based on the completed Hankel tensor. We present numerous simulation results in Section 5 and provide our conclusions in Section 6.

2. Tensor Signal Model with Array Sensor Failure

As shown in Figure 1, we consider an $M\times N$ URA

$$\mathbb{U}=\left\{(x_{\mathbb{U}},y_{\mathbb{U}})\right|x_{\mathbb{U}}=[0,M-1]d,y_{\mathbb{U}}=[0,N-1]d\},$$

(1)

where the unit inter-element spacing d equals to half of the signal wavelength. Assume that K narrowband far-field uncorrelated signals impinge on the URA $\mathbb{U}$ from directions $({\theta }_k,{\varphi }_k)$ for $k=1,2,\cdots ,K$, where ${\theta }_k$ and ${\varphi }_k$, respectively, denote the azimuth and elevation of the k-th source. For the conventional matrix-based signal processing methods, the array received signals at the t-th snapshot are modeled as a vector $\mathit{x}\left(t\right)\in {\mathbb{C}}^{MN}$, while the total T snapshots will then be stacked into a matrix

$$\mathit{X}=[\mathit{x}\left(1\right),\mathit{x}\left(2\right),\cdots ,\mathit{x}\left(T\right)]\in {\mathbb{C}}^{MN\times T}.$$

(2)

Here, $\mathbb{C}$ represents the set of complex numbers. Nevertheless, such operation fails to reflect the multi-spatial property of the URA signals, resulting in potential information loss.

Figure 1. The deployed URA with impaired sensors.

Figure 1. The deployed URA with impaired sensors.

To preserve the multi-spatial signal characteristics, we model the t-th snapshot as a matrix $\mathit{X}\left(t\right)\in {\mathbb{C}}^{M\times N}$, and then concatenate T snapshots to a 3D tensor, defined as

$$\mathcal{X}=\sum _{k=1}^K\mathit{a}\left({\mu }_k\right)\circ \mathit{a}\left({\nu }_k\right)\circ {\mathit{s}}_k+\mathcal{N}\in {\mathbb{C}}^{M\times N\times T}.$$

(3)

Here,

$$\begin{array}{cc}\hfill \mathit{a}\left({\mu }_k\right)& ={[1,e^{-j\pi {\mu }_k},\cdots ,e^{-j\pi (M-1){\mu }_k}]}^{\mathrm{T}},\hfill \\ \hfill \mathit{a}\left({\nu }_k\right)& ={[1,e^{-j\pi {\nu }_k},\cdots ,e^{-j\pi (N-1){\nu }_k}]}^{\mathrm{T}}\hfill \end{array}$$

(4)

are, respectively, the steering vectors of $\mathbb{U}$ along the x-axis and the y-axis with the angular projection factors

$$\begin{array}{cc}\hfill {\mu }_k& =sin{\varphi }_{k}cos{\theta }_k,\hfill \\ \hfill {\nu }_k& =sin{\varphi }_{k}sin{\theta }_k,\hfill \end{array}$$

(5)

${\mathit{s}}_k={[s_k\left(1\right),s_k\left(2\right),\cdots ,s_k\left(T\right)]}^{\mathrm{T}}$ is the signal waveform of the k-th source, ∘ denotes the outer product, $j=\sqrt{-1}$, ${(·)}^{\mathrm{T}}$ denotes the transpose operator, $\mathcal{N}$ is an additive Gaussian white noise tensor, i.e., its t-th slice $\mathit{N}\left(t\right)\sim \mathcal{C}\mathcal{N}(0,{\delta }_n^2\mathcal{I})$, $\forall t=1,2,\cdots ,T$, where ${\delta }_n^2$ denotes the noise power, and $\mathcal{I}$ denotes the identity tensor. In this way, the tensor signal model perfectly preserves the multi-spatial characteristics, which facilitates subsequent tensor optimizations for enhanced DOA estimation performance. Note that, the signal part of $\mathcal{X}$ in (3) follows a CP model, which is characterized by the sum of a certain number of rank-1 outer products.

The formulated tensor in (3) offers an ideal signal model for a completely functional URA. Unfortunately, when sensor failure occurs, the signals received at the impaired sensors will be missing. In this case, the corresponding URA signals can be characterized by an incomplete tensor

$$\tilde{\mathcal{X}}=\mathcal{M}⊛\mathcal{X}\in {\mathbb{C}}^{M\times N\times T},$$

(6)

where $\mathcal{M}\in {\{0,1\}}^{M\times N\times T}$ is the mask tensor indicating the positions of impaired sensors, and ⊛ denotes the Hadmard product. In particular, ${\mathcal{M}}_{(m,n,:)}=1$ if $(m,n)$ is the position of the unimpaired sensor, whereas ${\mathcal{M}}_{(m,n,:)}=0$ if $(m,n)$ is the position of the impaired sensor with $m=1,2,\cdots ,M$ and $n=1,2,\cdots ,N$.

As demonstrated in (6), each impaired sensor results in an entirely missing fiber in the corresponding position of the incomplete tensor $\tilde{\mathcal{X}}$. Moreover, as the number of impaired sensors increases, there might even exist whole missing slices in $\tilde{\mathcal{X}}$. As such, it is obvious that the sensor failure results in a structural-missing pattern for the tensor signal instead of the common random-missing pattern. The structural-missing tensor signal prevents a direct implementation of the conventional low-rank tensor completion, which pushes us to develop an effective structural-missing tensor completion method.

3. Structural-Missing Tensor Completion

In this section, we propose a structural-missing tensor completion method, facilitating the effective DOA estimation robust to sensor failures. The structural-missing tensor completion method includes spatio-temporal dimension augmentation and augmented tensor nuclear norm minimization optimization. The details in implementing the Hankel tensor formulation and solving the Hankel tensor completion problem are both provided.

3.1. Spatio-Temporal Dimension Augmentation for Missing Element Dispersion

Due to the existence of structurally missing elements in the incomplete tensor signal $\tilde{\mathcal{X}}$, the conventional low-rank regularization on an incomplete tensor with randomly missing elements becomes infeasible. Thus, to guarantee effective low-rank regularization on the tensor signal, we first design a spatio-temporal dimension augmentation method to disperse the structurally missing elements into the higher-dimensional space.

In particular, considering that different dimensions of the tensor signal present strong spatio-temporal relevance, properly augmenting these dimensions will enhance the capability of parameter retrieval, and also helps to disperse the structurally missing elements. The tensorial Hankelization approach is known for regularly embedding spatial/temporal delay information, which can be adopted for spatio-temporal dimension augmentation of the structural-missing tensor signal $\tilde{\mathcal{X}}$. To be more specific, as shown in Figure 2, $\tilde{\mathcal{X}}$ can be transformed into a six-dimensional (6D) dimensional-augmented Hankel tensor $\mathcal{H}\in {\mathbb{C}}^{{\tau }_1\times (M-{\tau }_1+1)\times {\tau }_2\times (N-{\tau }_2+1)\times {\tau }_3\times (T-{\tau }_3+1)}$, defined as

$$\mathcal{H}=fol{d}_{\left({\tau }_i\right)}\left(\tilde{\mathcal{X}}{\times }_1{\mathit{S}}_1{\times }_2{\mathit{S}}_2{\times }_3{\mathit{S}}_3\right),$$

(7)

where the folding operator $fol{d}_{\left({\tau }_i\right)}$ organizes a ${\tau }_1(L_1-{\tau }_1+1)\times {\tau }_2(L_2-{\tau }_2+2)\times \cdots \times {\tau }_I(L_I-{\tau }_I+1)$ tensor into another ${\tau }_1\times (L_1-{\tau }_1+1)\times {\tau }_2\times (L_2-{\tau }_2+1)\times \cdots \times {\tau }_I\times (L_I-{\tau }_I+1)$ tensor, ${\tau }_i$ denotes the duplication factor that expands the i-th dimension from the size of $L_i$ to ${\tau }_i(L_i-{\tau }_i+1)$, ${\times }_i$ denotes the mode-i tensor-matrix product, and ${\mathit{S}}_i\in {\left\{0,1\right\}}^{{\tau }_i(L_i-{\tau }_i+1)\times L_i}$ denotes the duplication matrices for $i=1,2,\cdots ,I$ with I being the tensor order (For the structural-missing tensor signal $\tilde{\mathcal{X}}$, we have $I=3$, but such an operation can be easily generalized to 3D cubic array with higher-dimensional tensor signals). Here, the $(p,q)$-th element of ${\mathit{S}}_i$ can be denoted by

$$\begin{array}{c}\hfill {\mathit{S}}_{i_{(p,q)}}=\left\{\begin{array}{cc}1\hfill & q=\lfloor p/{\tau }_i\rfloor +mod(p,{\tau }_i),\hfill \\ 0\hfill & \mathrm{others},\hfill \end{array}\right.\end{array}$$

(8)

where $\lfloor ·\rfloor$ denotes the rounding down operator, and $mod(·,·)$ denotes the residual operator. Since the size of the expanded dimensions should be a positive integer, we have ${\tau }_i<L_i$.

By performing spatio-temporal dimension augmentation of the structural-missing tensor signal as demonstrated in (7), its missing fibers or slices are now dispersed across six dimensions of the dimensional-augmented Hankel tensor $\mathcal{H}$, such that the low-rank tensor completion can be subsequently incorporated for its completion.

3.2. Dimensional-Augmented Hankel Tensor Completion

Thanks to the imposed spatio-temporal dimension augmentation, the resultant incomplete Hankel tensor $\mathcal{H}$ has sufficiently dispersed missing elements along with augmented spatio-temporal dimensions. As such, the augmented Hankel tensor $\mathcal{H}$ can be completed based on the low CP rank regularization (Note that, the Tucker-based tensor completion causes a problem of model mismatch in this scenario. Specifically, the Tucker decomposition yields a core tensor with uncertain Tucker ranks, which is inconsistent with the tensor signal model. Meanwhile, Tucker ranks should be manually initialized for Tucker decomposition, which causes performance deterioration if they are incorrectly set. Moreover, CPD based on the residual error minimization reconstructs a completed tensor, where the observed elements are modified. However, the observed elements in the tensor signal are expected to remain unchanged to avoid model deviations. As such, these tensor completion methods are not applicable to robust DOA estimation), which is characterized by the following optimization problem

$$\begin{array}{cc}\hfill & \underset{\overline{\mathcal{H}}}{min}{\mathrm{r}}_{\mathrm{CP}}\left(\overline{\mathcal{H}}\right)\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.{\mathrm{P}}_{\overline{\Omega }}\left(\overline{\mathcal{H}}\right)={\mathrm{P}}_{\overline{\Omega }}\left(\mathcal{H}\right),\hfill \end{array}$$

(9)

where the optimization variable $\overline{\mathcal{H}}\in {\mathbb{C}}^{{\tau }_1\times (M-{\tau }_1+1)\times {\tau }_2\times (N-{\tau }_2+1)\times {\tau }_3\times (T-{\tau }_3+1)}$ is the completed dimensional-augmented Hankel tensor, ${\mathrm{r}}_{\mathrm{CP}}$ denotes the CP rank of a tensor, $\Omega$ is the index set of missing elements in $\mathcal{H}$, and $\overline{\Omega }$ is the complement set of $\Omega$, i.e., the index set of the observed elements in $\mathcal{H}$. The proposed dimensional-augmented Hankel tensor completion problem in (9) minimizes the CP rank of $\overline{\mathcal{H}}$ under the constraint that all observed elements in $\mathcal{H}$ are kept in $\overline{\mathcal{H}}$. It is worth mentioning that the above optimization problem (9) is non-convex due to the NP-hardness of determining tensor ranks.

To solve the non-convex dimensional-augmented Hankel tensor completion problem (9), it is converted into a convex tensor nuclear norm minimization problem. The reason for this lies in that the nuclear norm is proved to be the tightest convex relaxation of the rank function. Specifically, the dimensional-augmented Hankel tensor nuclear norm minimization problem can be represented as

$$\begin{array}{cc}\hfill & \underset{\overline{\mathcal{H}}}{min}{\parallel \overline{\mathcal{H}}\parallel }_{*}=\sum _{j=1}^6{\alpha }_j{\parallel {\left[\overline{\mathcal{H}}\right]}_j\parallel }_{*}\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.{\mathrm{P}}_{\overline{\Omega }}\left(\overline{\mathcal{H}}\right)={\mathrm{P}}_{\overline{\Omega }}\left(\mathcal{H}\right),\hfill \end{array}$$

(10)

where ${\parallel ·\parallel }_{*}$ denotes the nuclear norm, ${[·]}_j$ denotes mode-j tensor unfolding, and ${\alpha }_j⩾0$ is the coefficient satisfying ${\sum }_{j=1}^6{\alpha }_j=1$. Note that, the matrix nuclear norm is proven to be the tightest relaxation of the rank function, and the tensor can be characterized by its unfolding matrices along respective dimensions. As such, the tensor nuclear norm serves as a convex combination of the nuclear norms of its matrix unfoldings, which makes it possible to solve the above problem with convex optimization approaches.

3.3. ADMM for Solving the Hankel Tensor Completion Problem

Normally, the optimization problem (10) is difficult to solve because the tensor nuclear norm terms ${\parallel \overline{\mathcal{H}}\parallel }_{*}$ are non-smooth. Hence, we adopt the alternating direction method of multipliers (ADMM), which can effectively solve large-scale optimization problems with non-smooth objectives [56]. The implementation of the ADMM for solving the Hankel tensor completion is elaborated in the following.

Specifically, in order to make sure that the six nuclear norm terms ${\parallel {\left[\overline{\mathcal{H}}\right]}_j\parallel }_{*}$ in the objective of the dimensional-augmented Hankel tensor nuclear norm minimization problem (10) can be independently optimized, we define six auxiliary tensors ${\mathcal{Y}}_j$ with identical initialization. As such, the dimensional-augmented Hankel tensor nuclear norm minimization problem (10) can be converted into an equivalent form

$$\begin{array}{cc}\hfill & \underset{\overline{\mathcal{H}}}{min}{\parallel \overline{\mathcal{H}}\parallel }_{*}=\sum _{j=1}^6{\alpha }_j{\parallel {\left[\overline{\mathcal{H}}\right]}_j\parallel }_{*}\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.{\mathrm{P}}_{\overline{\Omega }}\left(\overline{\mathcal{H}}\right)={\mathrm{P}}_{\overline{\Omega }}\left(\mathcal{H}\right),\hfill \\ \hfill & {\mathcal{Y}}_j-\overline{\mathcal{H}}=\mathcal{O},j=1,2,\cdots ,6,\hfill \end{array}$$

(11)

where $\mathcal{O}$ is zero tensor with appropriate size.

Then, we introduce six dual variables ${\mathcal{G}}_j\in {\mathbb{C}}^{{\tau }_1\times (M-{\tau }_1+1)\times {\tau }_2\times (N-{\tau }_2+1)\times {\tau }_3\times (T-{\tau }_3+1)}$ corresponding to the six dimensions of the optimization variables $\overline{\mathcal{H}}$ with $j=1,2,\cdots ,6$. Based on that, the augmented Lagrangian function of the dimensional-augmented Hankel tensor nuclear norm minimization problem (11) is defined as

$$\begin{array}{c}\hfill \mathcal{L}\left(\overline{\mathcal{H}},{\mathcal{Y}}_j,{\mathcal{G}}_j\right)={\parallel \overline{\mathcal{H}}\parallel }_{*}+<{\mathcal{Y}}_j-\overline{\mathcal{H}},{\mathcal{G}}_j>+\frac{\rho }{2}{∥{\mathcal{Y}}_j-\overline{\mathcal{H}}∥}_{\mathrm{F}}^2,\end{array}$$

(12)

where $\rho >0$ is the penalty constant, ${\parallel ·\parallel }_{\mathrm{F}}$ denotes the Frobenius norm, and $<·,·>$ denotes the inner product between two tensors. Subsequently, at the $(w+1)$-th iteration, the variables $\overline{\mathcal{H}}$, ${\mathcal{Y}}_j$ and ${\mathcal{G}}_j$ can be, respectively, updated as

$$\begin{array}{cc}\hfill {\mathcal{Y}}_j^{(w+1)}& =\arg \underset{{\mathcal{Y}}_j}{min}\mathcal{L}\left({\overline{\mathcal{H}}}^{\left(w\right)},{\mathcal{Y}}_j,{\mathcal{G}}_j^{\left(w\right)}\right),\hfill \\ \hfill {\overline{\mathcal{H}}}^{(w+1)}& =\arg \underset{\overline{\mathcal{H}}}{min}\mathcal{L}\left(\overline{\mathcal{H}},{\mathcal{Y}}_j^{(w+1)},{\mathcal{G}}_j^{\left(w\right)}\right),\hfill \\ \hfill {\mathcal{G}}_j^{(w+1)}& ={\mathcal{G}}_j^{\left(w\right)}-\rho \left({\mathcal{Y}}_j^{(w+1)}-{\overline{\mathcal{H}}}^{(w+1)}\right),\hfill \end{array}$$

(13)

where ${\overline{\mathcal{H}}}^{(w+1)}$ represents the optimization variable at the $(w+1)$-th iteration.

Using the principle of gradient descent to solve (13), we can obtain the closed-form solution to ${\mathcal{Y}}_j^{(w+1)}$ and ${\overline{\mathcal{H}}}^{(w+1)}$, i.e.,

$$\begin{array}{cc}\hfill & {\mathcal{Y}}_j^{(w+1)}={\mathrm{fold}}_{\left(j\right)}\left({{\rm Y}}_{\frac{{\alpha }_j}{\rho }}\left({\left[{\overline{\mathcal{H}}}^{\left(w\right)}\right]}_j+\frac{1}{\rho }{\left[{\mathcal{G}}_j^{(w+1)}\right]}_j\right)\right),\hfill \\ \hfill & {\mathrm{P}}_{\Omega }\left({\overline{\mathcal{H}}}^{(w+1)}\right)={\mathrm{P}}_{\Omega }\left(\frac{1}{6}\left(\sum _{j=1}^6{\mathcal{Y}}_j^{(w+1)}-\frac{1}{\rho }{\mathcal{G}}_j^{\left(w\right)}\right)\right),\hfill \\ \hfill & {\mathrm{P}}_{\overline{\Omega }}\left({\overline{\mathcal{H}}}^{(w+1)}\right)={\mathrm{P}}_{\overline{\Omega }}\left(\mathcal{H}\right),\hfill \end{array}$$

(14)

where

$${{\rm Y}}_{\frac{{\alpha }_j}{\rho }}\left(\mathit{U}\right)={\mathit{F}}_U{\Sigma }_U^{\left(\frac{{\alpha }_j}{\rho }\right)}{\mathit{W}}_U$$

(15)

denotes the shrinkage singular value decomposition of a matrix $\mathit{U}\in {\mathbb{C}}^{U_1\times U_2}$, ${\mathit{F}}_U$, ${\mathit{W}}_U$ represent the left singular matrix and the right singular matrix, respectively, and

$${\Sigma }_U^{\left(\frac{{\alpha }_j}{\rho }\right)}=\mathrm{diag}\left(max\left({\eta }_l-\frac{{\alpha }_j}{\rho },0\right)\right).$$

(16)

Here, ${\eta }_l$ denotes the singular value of $\mathit{U}$, $l\in \left\{1,2,\cdots ,min(U_1,U_2)\right\}$, $min(·)$ and $max(·)$, respectively, denote the minimum and maximum value, and $\mathrm{diag}(·)$ forms a diagonal matrix from its arguments.

By repeating the above updating procedure, we can finally obtain the optimized variable, i.e., the completed Hankel tensor $\overline{\mathcal{H}}$. In particular, the ADMM for solving the Hankel tensor completion problem converges when the relative error of the completed Hankel tensor between two iterations is less than a threshold $\xi >0$, that is,

$$\frac{{\parallel {\overline{\mathcal{H}}}^{(w+1)}-{\overline{\mathcal{H}}}^{\left(w\right)}\parallel }_{\mathrm{F}}}{{\parallel {\overline{\mathcal{H}}}^{\left(w\right)}\parallel }_{\mathrm{F}}}⩽\xi .$$

(17)

Note that, the ADMM for solving the general tensor nuclear norm optimization problem is proven to be globally convergent. The detailed proof is provided in [43].

4. Robust DOA Estimation with Array Sensor Failure

In this section, we propose an inverse Hankelization approach to derive the completed covariance tensor based on the designed completed dimensional-augmented Hankel tensor. The covariance tensor is then decomposed to obtain estimated DOAs with high robustness to sensor failure. We also present analysis on the computational complexity for the proposed algorithm.

4.1. Inverse Hankelization for Completed Covariance Tensor Derivation

Based on the completed 6D dimensional-augmented Hankel tensor $\overline{\mathcal{H}}$, we propose the inverse Hankelization to recover a completed 3D tensor signal $\mathcal{Y}$ which serves as the equivalent signals received at an unimpaired URA $\mathbb{U}$. The reduced dimensionality not only contributes to higher efficiency, but also facilitates the calculation of tensor statistics. Based on that, the completed covariance tensor of $\mathcal{Y}$ can be derived for DOA estimation.

As shown in Figure 2, according to the spatio-temporal dimension augmentation of the incomplete tensor signal which involves dimensional duplication and folding, its inverse procedure can be represented as

$$\mathcal{Y}=unfol{d}_{\left({\tau }_i\right)}\left(\overline{\mathcal{H}}\right){\times }_1{\mathit{S}}_1^{†}{\times }_2{\mathit{S}}_2^{†}{\times }_3{\mathit{S}}_3^{†}\in {\mathbb{C}}^{M\times N\times T},$$

(18)

where ${(·)}^{†}$ represents the Moore–Penrose pseudo-inverse operator, and

$$unfol{d}_{\left({\tau }_i\right)}(·)=fol{d}_{\left({\tau }_i\right)}^{-1}(·)$$

(19)

represents the inverse operator of $fol{d}_{\left({\tau }_i\right)}(·)$. Note that, the completed tensor signal $\mathcal{Y}$ is an approximate of the ideal tensor signal $\mathcal{X}$ in (3). Thus, $\mathcal{Y}$ should also conform to the CP representation as $\mathcal{X}$.

Furthermore, the second-order tensor statistics of the URA $\mathbb{U}$ can be derived from the completed tensor signal $\mathcal{Y}$ for DOA estimation. Specifically, we calculate the four-dimensional (4D) covariance tensor of $\mathcal{Y}$ with a CP representation as

$$\begin{array}{cc}\hfill \mathcal{R}& =\mathrm{E}\left\{\frac{1}{T}{〈\mathcal{Y},{\mathcal{Y}}^{*}〉}_3\right\}\in {\mathbb{C}}^{M\times N\times M\times N}\hfill \\ \hfill & =\sum _{k=1}^K{\sigma }_{s_k}^2\mathit{a}\left({\mu }_k\right)\circ \mathit{a}\left({\nu }_k\right)\circ {\mathit{a}}^{*}\left({\mu }_k\right)\circ {\mathit{a}}^{*}\left({\nu }_k\right)+{\sigma }_n^2\mathcal{I},\hfill \end{array}$$

(20)

where ${\sigma }_{s_k}^2=\mathrm{E}\left\{s_k\left(t\right)s_k^{*}\left(t\right)\right\}$ denotes the power of k-th source signal, ${\sigma }_n^2$ is the noise power, $\mathcal{I}$ denotes the 4D identity tensor, ${(·)}^{*}$ denotes the conjugation operator, $\mathrm{E}\left\{·\right\}$ denotes the statistical expectation, and ${〈·,·〉}_i$ denotes the tensor contraction along the i-th dimension. Now, $\mathcal{R}$ can be regarded as the completed covariance tensor of the unimpaired URA $\mathbb{U}$. Thus, we can decompose it to retrieve the angular parameters embedded in the CP factors of $\mathcal{R}$, leading to effective DOA estimation under sensor failure. The implementation of decomposition $\mathcal{R}$ is elaborated in the next subsection. In practice, the covariance tensor can be approximated by its sample version

$$\widehat{\mathcal{R}}=\frac{1}{T}{〈\mathcal{Y},{\mathcal{Y}}^{*}〉}_3$$

(21)

due to the finite number of snapshots.

4.2. Covariance Tensor CPD for Robust DOA Estimation

By decomposing the covariance tensor $\mathcal{R}$ to retrieve the estimated steering vectors $\mathit{a}\left({\widehat{\mu }}_k\right)$ and $\mathit{a}\left({\widehat{\nu }}_k\right)$, the estimated angular projection factors ${\widehat{\mu }}_k$ and ${\widehat{\nu }}_k$ can be obtained, leading to the closed-form solution of both azimuth and elevation. Considering that $\mathcal{R}$ is represented in a CP model, we apply the CPD to $\mathcal{R}$ by solving the following least squares problem

$$\underset{\widehat{\mathit{A}}\left(\mu \right),\widehat{\mathit{A}}\left(\nu \right)}{min}{∥\mathcal{R}-\mathcal{J}{\times }_1\widehat{\mathit{A}}\left(\mu \right){\times }_2\widehat{\mathit{A}}\left(\nu \right){\times }_3{\widehat{\mathit{A}}}^{*}\left(\mu \right){\times }_4{\widehat{\mathit{A}}}^{*}\left(\nu \right)∥}_{\mathrm{F}}^2,$$

(22)

where

$$\begin{array}{cc}\hfill \widehat{\mathit{A}}\left(\mu \right)& =\left[\mathit{a}\left({\widehat{\mu }}_1\right),\mathit{a}\left({\widehat{\mu }}_2\right),\cdots ,\mathit{a}\left({\widehat{\mu }}_K\right)\right]\in {\mathbb{C}}^{M\times K},\hfill \\ \hfill \widehat{\mathit{A}}\left(\nu \right)& =\left[\mathit{a}\left({\widehat{\nu }}_1\right),\mathit{a}\left({\widehat{\nu }}_2\right),\cdots ,\mathit{a}\left({\widehat{\nu }}_K\right)\right]\in {\mathbb{C}}^{N\times K},\hfill \end{array}$$

(23)

are the estimated steering matrices of $\mathbb{U}$, $\mathcal{J}\in {\mathbb{C}}^{K\times K\times K}$ is the signal power tensor with ${\sigma }_{s_k}^2$ on its main diagonal.

The above least squares problem (22) can be solved by iteratively updating $\widehat{\mathit{A}}\left(\mu \right)$ and $\widehat{\mathit{A}}\left(\nu \right)$ as

$$\begin{array}{cc}\hfill \widehat{\mathit{A}}\left(\mu \right)& =\arg \underset{\widehat{\mathit{A}}\left(\mu \right)}{min}{∥{\left[\mathcal{R}\right]}_1-\widehat{\mathit{A}}\left(\mu \right){\left({\widehat{\mathit{A}}}^{*}\left(\nu \right)\odot {\widehat{\mathit{A}}}^{*}\left(\mu \right)\odot \widehat{\mathit{A}}\left(\nu \right)\right)}^{\mathrm{T}}∥}_{\mathrm{F}}^2,\hfill \\ \hfill \widehat{\mathit{A}}\left(\nu \right)& =\arg \underset{\widehat{\mathit{A}}\left(\nu \right)}{min}{∥{\left[\mathcal{R}\right]}_2-\widehat{\mathit{A}}\left(\nu \right){\left({\widehat{\mathit{A}}}^{*}\left(\nu \right)\odot {\widehat{\mathit{A}}}^{*}\left(\mu \right)\odot \widehat{\mathit{A}}\left(\mu \right)\right)}^{\mathrm{T}}∥}_{\mathrm{F}}^2,\hfill \end{array}$$

(24)

where ⊙ denotes the Khatri–Rao product. For each iteration, the closed-form solution of the estimated steering matrices can be computed as

$$\begin{array}{cc}\hfill \widehat{\mathit{A}}\left(\mu \right)& ={\left[\mathcal{R}\right]}_1{\left[{\left({\widehat{\mathit{A}}}^{*}\left(\nu \right)\odot {\widehat{\mathit{A}}}^{*}\left(\mu \right)\odot \widehat{\mathit{A}}\left(\nu \right)\right)}^{\mathrm{T}}\right]}^{†},\hfill \\ \hfill \widehat{\mathit{A}}\left(\nu \right)& ={\left[\mathcal{R}\right]}_2{\left[{\left({\widehat{\mathit{A}}}^{*}\left(\nu \right)\odot {\widehat{\mathit{A}}}^{*}\left(\mu \right)\odot \widehat{\mathit{A}}\left(\mu \right)\right)}^{\mathrm{T}}\right]}^{†}.\hfill \end{array}$$

(25)

The iterations repeat until the relative error of the decomposed covariance tensor between successive iterations is smaller than a convergence threshold, and the resulting estimation of steering matrices can be used for retrieving the angular parameters.

As such, the estimation of the angular projection factors ${\widehat{\mu }}_k$ and ${\widehat{\nu }}_k$ can be obtained as

$$\begin{array}{cc}\hfill {\widehat{\mu }}_k& =\frac{1}{M}\sum _{m=1}^M\angle ({\mathit{a}}_{m+1}\left({\widehat{\mu }}_k\right)/{\mathit{a}}_m\left({\widehat{\mu }}_k\right))/\pi ,\hfill \\ \hfill {\widehat{\nu }}_k& =\frac{1}{N}\sum _{n=1}^N\angle ({\mathit{a}}_{n+1}\left({\widehat{\nu }}_k\right)/{\mathit{a}}_n\left({\widehat{\nu }}_k\right))/\pi ,\hfill \end{array}$$

(26)

where ${\mathit{a}}_m\left({\widehat{\mu }}_k\right)$, ${\mathit{a}}_n\left({\widehat{\nu }}_k\right)$, respectively, represent the m-th element and n-th element of $\widehat{\mathit{a}}\left({\mu }_k\right)$, $\widehat{\mathit{a}}\left({\nu }_k\right)$, and ∠ denotes the phase of a complex number. Finally, according to the relationship between angular projection factors $({\mu }_k,{\nu }_k)$ and $({\theta }_k,{\varphi }_k)$ as defined in (5), the azimuth and elevation of the k-th source can be calculated in closed-forms as

$$\begin{array}{cc}\hfill {\widehat{\theta }}_k& =\arctan ({\widehat{\nu }}_k/{\widehat{\mu }}_k),\hfill \\ \hfill {\widehat{\varphi }}_k& =\arcsin \left(\sqrt{{\widehat{\nu }}_k^2+{\widehat{\mu }}_k^2}\right).\hfill \end{array}$$

(27)

4.3. Analysis on Computational Complexity

In this subsection, we provide theoretical analysis on the computational complexity of the proposed algorithm. The proposed algorithm mainly involves the Hankel tensor completion and the completed covariance tensor decomposition, whose computational complexities are $\mathcal{O}\left((M+N+T)(M-{\tau }_1)(N-{\tau }_2)(T-{\tau }_3){\tau }_1{\tau }_2{\tau }_3{Q}_{\mathrm{TC}}\right)$ and $\mathcal{O}\left(M^2{N}^{2}T+(K^3+M^2{N}^2)K{Q}_{\mathrm{CPD}}\right)$. Here, $Q_{\mathrm{TC}}$ and $Q_{\mathrm{CPD}}$ are, respectively, the numbers of iterations for implementing the tensor completion and CPD. The computational complexity of the proposed algorithm can then be measured by $\mathcal{O}\left((M+N+T)(M-{\tau }_1)(N-{\tau }_2)(T-{\tau }_3){\tau }_1{\tau }_2{\tau }_3{Q}_{\mathrm{TC}}+M^2{N}^2(Q_{\mathrm{CPD}}+T)\right)$. Due to the globally convergent ADMM for tensor completion, the proposed Hankel tensor completion approach presents satisfactory convergence behavior. It normally converges within dozens of iterations, making $Q_{\mathrm{TC}}$ relatively small.

The computational complexity of the matrix completion (MC) method [23] is $\mathcal{Q}\left(T+log\left(\frac{1}{\eta }\right)(rMNT+r^2(MN+T))+T(M^2{N}^2+MN)+MN{Q}_{\mathrm{SG}}\right)$, where $\eta$ and $Q_{\mathrm{SG}}$, respectively, denote the convergence threshold and number of spectral searching grids for implementing the matrix completion. In addition, the computational complexity of the direct CPD method [37] is $\mathcal{O}\left(M^2{N}^{2}T+(K^3+M^2{N}^2)K{Q}_{\mathrm{CPD}}\right)$. The proposed algorithm is more computational costly compared to the direct CPD method due to the Hankel tensor completion. Also, the high-order tensor formulation and multi-linear tensor completion would cause increased complexities than the matrix completion and spectral searching procedures. However, the satisfactory convergence of the ADMM and closed-form solution of CPD enables comparable efficiency. Thus, though incorporating the high-order tensor formulation and optimization, the computational complexity of this proposed algorithm can be moderate compared with existing methods.

5. Simulation

In the simulation, we set the duplication factors $\{{\tau }_i,i=1,2,3\}$ for spatio-temporal dimension augmentation to $[3,3,10]$, and set the coefficients $\{{\alpha }_j,j=1,2,\cdots ,6\}$ for defining the Hankel tensor nuclear norm equally to $\frac{1}{6}$. We set the penalty constant of the ADMM augmented Lagrangian function (12) to $\rho ={10}^{-3}$, and set the threshold of ADMM to $\xi ={10}^{-12}$. Since the positions of impaired sensors in the URA are unpredictable in real applications, they are set randomly for each simulated scenario. The root mean square error (RMSE)

$$\mathrm{RMSE}=\sqrt{\frac{1}{2K{N}_{\mathrm{MC}}}\sum _{n_{\mathrm{MC}}=1}^{N_{\mathrm{MC}}}\sum _{k=1}^K\left[{({\widehat{\theta }}_{k,n_{\mathrm{MC}}}-{\theta }_k)}^2+{({\widehat{\varphi }}_{k,n_{\mathrm{MC}}}-{\varphi }_k)}^2\right]}$$

(28)

is adopted as the evaluation metric, where $({\widehat{\theta }}_{k,n_{\mathrm{MC}}},{\widehat{\varphi }}_{k,n_{\mathrm{MC}}})$ is the estimate of $({\theta }_k,{\varphi }_k)$ for the $n_{\mathrm{MC}}$-th Monte Carlo trial, and $N_{\mathrm{MC}}=1000$ Monte Carlo trials are run for each data point to obtain the RMSE curves. Unless otherwise specified, we deploy a URA $\mathbb{U}$ with $M=7$ and $N=7$, and assume $K=2$ sources from the directions $({\theta }_1,{\varphi }_1)=({10}^{\circ },{10}^{\circ })$ and $({\theta }_2,{\varphi }_2)=({25}^{\circ },{25}^{\circ })$, respectively. The proposed algorithm is implemented with Tensorlab 3.0 on MATLAB R2021b, which develops useful tools for implementing the tensor formulation and CPD [57]. Meanwhile, the ADMM for solving the Hankel tensor completion problem can be implemented with well-studied tool in the literature [43].

5.1. Performance of the Proposed Algorithm with Different Array Geometries

We first evaluate the performance of the proposed algorithm with different array geometries. As claimed before, the $7\times 7$ URA is used as the benchmark. Then, we deploy two larger URAs with the sizes of $10\times 10$ and $12\times 12$. $K=2$ sources come from $({35}^{\circ },{25}^{\circ })$ and $({68}^{\circ },{57}^{\circ })$, respectively. The numbers of impaired sensors in these URAs are all set to 5. The signal-to-noise ratios (SNR) of sources are equally set to 10 dB, and the number of snapshots is fixed to $T=100$. We run 10 trials, and present the estimation results in Figure 3. It is clear that all URAs with different geometries can accurately estimate the source DOAs under sensor failure. Meanwhile, with a larger array size, the estimation accuracy increases.

5.2. Comparison of the Estimation Accuracy with Different Number of Damaged Sensors

In Figure 4, we compare the estimation accuracy of the proposed structural-missing tensor completion algorithm with those of the matrix completion method [23] and the CPD method [37], where the number of impaired sensors varies from 1 to 23. The SNRs of sources are equally set to 15 dB, and the number of snapshots is fixed to $T=300$.

It is clear that the proposed algorithm is superior to the compared methods in all simulated sensor failure scenarios. The improvement comes from the utilization of multi-spatial signal characteristics, and more importantly, the effective recovery of structural-missing elements in the formulated tensor signal. It is also worth mentioning that as the number of impaired sensors increases (30% sensor failure with 15 impaired sensors), the performance of all methods degrade, while the proposed algorithm presents a moderate reduction in estimation accuracy. The proposed method fails to estimate the DOA when the number of impaired sensors reaches 23 (47% of the sensors are impaired), which can be regarded as the boundary for the proposed algorithm. However, in practice, it is unusual to have a very high percentage of impaired sensors. The above findings demonstrate the high robustness of the proposed algorithm against serious sensor failure conditions.

Figure 4. Comparison of accuracy with different number of impaired sensors.

Figure 4. Comparison of accuracy with different number of impaired sensors.

5.3. Comparison of the Estimation Accuracy with Different SNRs and Number of Snapshots

Furthermore, the estimation accuracy of the above-mentioned methods is compared with 5 impaired sensors in Figure 5. The SNR varies whereas the number of snapshots is set to $T=100$ in Figure 5a, while the number of snapshots varies whereas the SNRs of sources are set to 15 dB in Figure 5b. The implementation of CPD on an unimpaired URA is presented as reference. In addition, considering the CRB is a common lower bound on the performance of DOA estimators, we use it as the benchmark in the simulation. Specifically, we calculate the weighted sum of CRB for a URA with randomly impaired sensors [30].

As shown in Figure 5a, the estimation accuracy of the proposed structural-missing tensor completion algorithm is the highest among the entire tested SNR regime as compared to other methods with array sensor failure. The proposed algorithm outperforms the CPD method with array sensor failure, because it effectively completes the tensor signals with missing fibers or slices. In contrast, the direct CPD on the incomplete tensor statistics suffers from worse capability of parameter retrieval. On the other hand, since the proposed tensor signal model preserves the original 2D angular information, the proposed algorithm also enjoys a significant improvement in accuracy compared to the conventional matrix completion method. Moreover, as the SNR increases, the RMSE curve of the proposed algorithm gradually approaches that of the CPD method for the unimpaired URA, and the two RMSE curves become very close with SNR higher than 10 dB, which also verifies the advantage of the proposed algorithm in overcoming the effect of impaired sensors. Similar result can be found in Figure 5b, where the estimation accuracy of the proposed algorithm is very close to that of the CPD method for the unimpaired URA, and is higher than all other methods.

5.4. Comparison of the Angular Resolution

In Figure 6, we compare the angular resolution of the evaluated methods with 5 impaired sensors, where the two sources $({\theta }_1,{\varphi }_1)$ and $({\theta }_2,{\varphi }_2)$ are placed closely. Here, the source $({\theta }_1,{\varphi }_1)$ is fixed at $({20}^{\circ },{40}^{\circ })$, and $({\theta }_2,{\varphi }_2)$ is generated with an angular spacing of

$$\delta =\sqrt{{|{\theta }_2-{\theta }_1|}^2+{|{\varphi }_2-{\varphi }_1|}^2}$$

(29)

to $({\theta }_1,{\varphi }_1)$. The evaluated methods are identify to distinguish the two sources if

$$\sqrt{{|{\widehat{\theta }}_{k,n_{\mathrm{MC}}}-{\theta }_k|}^2+{|{\widehat{\varphi }}_{k,n_{\mathrm{MC}}}-{\varphi }_k|}^2}<\delta$$

(30)

for each trial. The successful rate will be calculated from the number of successful trials over the $N_{\mathrm{MC}}=1000$ Monte Carlo trials. The SNRs of sources are set to 15 dB, the number of snapshots is fixed at $T=100$, and the angular spacing varies from $2^{\circ }$ to $8^{\circ }$.

Obviously, the angular resolution of the proposed algorithm is higher than the other methods. Such improvement is benefited from the enforced structural-missing tensor completion, which enables to utilize the entire URA aperture for high-resolution DOA estimation. In contrast, for other methods, the URA aperture is not fully exploited due to the impaired sensors, resulting in degraded angular resolution of DOA estimation. Note that the proposed algorithm can still identify source DOAs with the angular spacing of only $2^{\circ }$, whereas the compared methods almost fail.

5.5. Performance of Source Identifiability

The source identifiability of using the 4D covariance tensor $\mathcal{R}\in {\mathbb{C}}^{M\times N\times M\times N}$ corresponding to the $M\times N$ URA $\mathbb{U}$ can be measured by the Kruskal’s condition. Generally speaking, for the 4D tensor $\mathcal{A}\in {\mathbb{C}}^{I_1\times I_2\times I_3\times I_4}$ with rank R, its CPD is essentially unique if

$$\begin{array}{c}\hfill R\le I_4,R(R-1)\le I_1{I}_2{I}_3(3I_1{I}_2{I}_3-I_1{I}_2-I_2{I}_3-I_1{I}_3-I_1-I_2-I_3+3)/4.\end{array}$$

(31)

As for the 4D covariance tensor $\mathcal{R}$, the maximum of R can be calculated as 7. This means that a $7\times 7$ URA can identify up to $K=7$ sources.

Based on the above analysis, we validate the source identifiability of the proposed algorithm in Figure 7, where the azimuth and elevation angles of 7 sources are uniformly distributed within $[-{50}^{\circ },{70}^{\circ }]$. The SNRs of sources are set to 15 dB and the number of snapshots is fixed at $T=200$. The proposed algorithm can locate all these sources with high accuracy. Such a finding verifies the source identifiability of the proposed algorithm.

In addition, we test the performance of the proposed algorithm for 8 sources in Figure 8 with the same simulation settings. Although the analytical result presented in (31) offers the theoretical bound on the maximum number of distinguishable sources, the proposed algorithm is still capable of estimating more sources with slightly lower accuracy, which further highlights its robustness in practical usage.

Figure 7. Source identifiability of the proposed algorithm.

Figure 7. Source identifiability of the proposed algorithm.

Figure 8. Estimation of 8 sources.

Figure 8. Estimation of 8 sources.

5.6. Comparison with Other Completion Methods

To validate the effectiveness of the proposed structural-missing completion algorithm with serious sensor failure, in Figure 9, we compare the estimation accuracy of the proposed algorithm with those of the CPD-based completion method and Tucker-based completion method, where the number of impaired sensors reaches 10, i.e., 20% of the sensors are impaired. As illustrated in the Introduction, the CPD-based completion method use the CP model to approximate the completed tensor, whereas the observed elements will also be modified. The Tucker-based completion method imposes the Tucker rank increment to complete the missing fibers, but requires exhausting rank searching and does not fit the tensor signal model.

It is clear that the estimation accuracy of the proposed algorithm is superior to the other two completion methods. It is worth noting that the error between the manually set Tucker rank and the real Tucker rank will be amplified as the number of failed sensors increases, which results in degraded performance for the Tucker-based completion method. This highlights the advantage of the enforced low CP rank regularization for structural-missing tensor signal completion.

5.7. Comparison of the Computational Time

Finally, in Figure 10, we compare the computational time of the proposed algorithm with those of the matrix completion method and the direct CPD method. The simulation settings are kept the same as in the former simulations. As demonstrated in Section 4.3, the complexity of the proposed algorithm is greater than those of the direct CPD method and the matrix completion method, due to the high-order tensor formulation and the structural-missing tensor completion for performance enhancement. Therefore, the proposed algorithm requires more computational time in pursuit of higher estimation accuracy and resolution as shown in the previous simulation. However, thanks to the efficient convergence of tensor completion and closed-form tensor decomposition, the proposed algorithm is capable of maintaining moderate efficiency. It is clear that the computational time of the proposed algorithm and the matrix completion method is at a similar level. Therefore, the proposed algorithm can still offer practical DOA estimation with significantly enhanced performance.

6. Conclusions

In this paper, a structural-missing tensor completion-based robust DOA estimation algorithm is proposed to address the array sensor failure problem. The structural-missing pattern of the incomplete tensor signal is investigated, based on which the spatio-temporal dimension augmentation method is designed to construct a six-dimensional augmented Hankel tensor with sufficiently dispersed missing elements. As such, the augmented Hankel tensor can be effectively completed by solving the Hankel tensor nuclear norm minimization problem with ADMM. Implementing inverse Hankelization on the completed Hankel tensor yields equivalent tensor signal of an unimpaired sensor array, such that the completed covariance tensor can be derived. Finally, the CPD of the covariance tensor leads to closed-form DOA estimation. According to both analytical and numerical results, the proposed algorithm achieves superior estimation performance compared to the competing methods, including enhanced estimation accuracy, angular resolution, robustness to sensor failure, source identifiability, and moderate computational complexity.