Direction of Arrival Estimation of Acoustic Sources with Unmanned Underwater Vehicle Swarm via Matrix Completion

With the rapid development of vibration and noise reduction technologies, underwater target detection is facing great challenges. Particularly, the task of high-resolution direction of arrival (DOA) estimation with sonar array is becoming more and more tough. In recent years, unmanned underwater vehicles (UUVs) have been developed considerably, with the improvements of target localization performance in terms of adaptability, detection range, operation efficiency, and antiinterference ability. Nevertheless, in general, the size of UUV is small such that current passive sonar systems usually have relatively limited localization accuracy, detection distance, and environmental robustness in complex ocean noise. This motivates us to present a new approach to construct a large-aperture virtual array with multiple small-aperture arrays of unmanned underwater vehicle swarm (UUVS) which consists of multiple UUVs in this paper. However, for the UUVS array, the received data could suffer from unobserved and corrupted samples. This makes it challenging to analyze and process large-aperture array data. Towards this end, the matrix completion technique is employed to recover the unobserved and corrupted data for virtual array construction based on the low rank property of array data matrix. The recovered matrix is then exploited for underwater target bearing estimation using the traditional DOA estimation approach. Numerical results verify that the proposed method is capable of detecting underwater targets with high precision and resolution.


Direction of Arrival Estimation of Acoustic Sources
As the basis of target localization, direction of arrival (DOA) estimation is an important research topic in underwater acoustic engineering and a prerequisite for antisubmarine, torpedo defense, and underwater acoustic countermeasure. During the past several decades, numerous DOA estimation methods have been developed and applied to acoustic sources [1]. Representative methods include subspace-based approaches like multiple signal classification (MUSIC), estimation of parameter via rotational invariance technique (ESPRIT), and their variants [2], as well as beamforming approaches such as the conventional beamforming (CBF), minimum variance distortionless response (MVDR), delay-and-sum, and their modifications [3].
Facing the increasingly complex ocean environment and various interference sources deployed, enhanced DOA estimation methods are desired to overcome the problem of sound fluctuations caused by the space-time nonuniform propagation of sound field [4] for precise target localization. For instance, in [5], a high-precision DOA estimation method for underwater acoustic sources is developed by exploiting the signal self-cancellation phenomenon of the MVDR beamformer. In [6], a maximum-likelihood (ML) DOA estimator for mixed signals (containing both sinusoidal and random components) from ships, submarines, or torpedoes is presented and shown to be superior to typical stochastic estimator assuming zero-mean Gaussian signals. More recently, a passive DOA estimation algorithm which is able to deal with closely-spaced multipath underwater signals is devised based on spatial time-frequency distributions (STFDs) in [7].

Underwater Unmanned Vehicle Swarm (UUVS)
To perform DOA estimation of underwater acoustic signals and other missions, the acoustic sensor arrays are now usually carried by unmanned underwater vehicles (UUVs) [8], which have recently been developed and adopted at an increasing rate. UUVs use remote control or automatic control to navigate in the underwater environment. They mainly refer to smart equipment that can replace divers or manned mini-submarine for detection, life-saving, torpedo removal, and other high-risk underwater operations. UUVs are also called robotic submarines or underwater robots, which can be mainly categorized as remote operated vehicles (ROVs) and autonomous underwater vehicles (AUVs) [8][9][10][11][12]. They normally use the submarine or surface ship as support platforms, and are equipped with a variety of sensors, energy and navigation equipment, and special weapon, etc.
UUVs could conduct electronic attacks in littorals by taking advantage of their small size and stealth to operate in areas unreachable by other platforms [13]. Nevertheless, it is known that the limited size of an UUV could impose restrictions on the array aperture and hence the detection performance. Thus, marine equipment is being developed towards the swarm intelligence. Compared with the traditional single-robot system, a swarm of UUVs is able to provide richer spatial-temporal information, which can be used for underwater mapping, exploration, and target tracking [14].
Swarming behavior, which originated in the biological world, is a collective behavior of organisms, such as insects, birds, fish, humans, and bacteria. Swarming technology of robotic is conceptually similar to fish swarm, bird swarm, bee swarm, and ant swarm. For the development of electronic information in the future, a variety of unmanned equipment will be applied to form an unmanned swarm system in the ocean [15], which can meet a variety of complex tasks. The marine equipment of unmanned swarm systems mainly includes unmanned aerial vehicle swarm (UAVS) [16], unmanned surface vessel swarm (USVS) [17], and unmanned underwater vehicle swarm (UUVS) [18]. In this paper, we focus on UUVS, and the schematic diagram is shown in Figure 1. Note that UUVs in the swarm can move with respect to each other in practice. Herein, we only consider a formation without internal movements.

Challenges of UUVS DOA Estimation
As mentioned above, compared with a single UUV, the UUVS could obtain better detection performance through the cooperation of multiple UUVs, and have various advantages such as strong adaptability, large detection range, high operation efficiency, and strong anti-interference ability [19]. The swarming capability of unmanned systems has been demonstrated by the SwarmDiver, manufactured by Aquabotix Technology Corporation, which can dive as deep as 600 feet while working in synch with other UUVs submersibles [20]. Note that, at present, the research on UUVS mainly focuses on network communication, formation, navigation, and the concept exploration of cooperative combat. Due to the complexity of the underwater acoustic field, the application of UUVS to acoustic passive DOA estimation is an emerging research topic. Nevertheless, there still exist many challenges in this area.
According to array signal processing theories, the resolution capability of DOA estimation approaches is dependent on the array aperture and number of sensors. More specifically, a larger aperture leads to higher resolution [21]. Thus, large apertures are preferred to achieve better estimation performance. However, in practical implementations, enlarging the array aperture may be difficult or even impossible for UUVs which have limited size to ensure their mobility and flexibility. Even the UUVS with multiple short linear arrays (as illustrated in Figure 2) can form a large-aperture array, the number of sensors may be still limited, and the resulting array is sparse. Accordingly, the DOA estimation performance could be unsatisfactory, mainly manifested as wide mainlobe width, low estimation accuracy, and resolution, etc. Various approaches have been proposed to expand array aperture, such as the synthetic aperture method [22] and a high-order based method [23]. However, these methods are generally not straightforwardly applicable to UUVS.

UUV #1
UUV #2 UUV #3 Far-field Source Another challenge for DOA estimation with UUVS array is that the received data are susceptible to missing or corrupted samples. For an UUVS, each UUV node collects measurements and uploads the data to a central place for further processing including compilation, modeling, and analysis. In practice, it may be impossible to obtain every measurement from every point in the network due to resource constraints, node outages, etc. [24]. Moreover, data missing could be caused by sensor failure [25]. However, it is usually difficult to carry out underwater equipment maintenance. The system performance may be significantly degraded if these aspects are ignored. From the signal processing point of view, the problem of sensor failure can be addressed by developing new algorithms that are functional in the presence of sensor failure and analyzing the robustness against sensor failure or data missing [26]. However, the performance could be limited.

Contributions of This Work
Motivated by the aforementioned challenges, in this work, we propose to use matrix completion (MC) techniques to deal with the geometry sparsity of the UUVS array and the issues of unobserved/corrupted entries. Specifically, we use MC to recover the unobserved data of the virtual sensors between two UUVs, such that a large-aperture virtual array can be obtained. Since the data matrix with unobserved columns cannot be directly recovered by MC, a method of structured MC with Hankel matrix transformation is presented. Furthermore, to deal with corrupted entries, a sparse matrix is introduced to the data matrix to take the data pollution into account. Accordingly, a robust matrix completion (RMC) algorithm is given to solve the resulting problem. With the reconstructed UUVS array data matrix, DOA estimation can be performed with traditional approaches.

UUVS Array Signal Model
As illustrated in Figure 2, a UUVS composed of U UUVs is considered, and each UUV is equipped with a P-element hydrophone uniform linear array with inter-element spacing d, then the number of elements in the receiving array (composed of U subarrays) is N 0 = UP. Suppose that free-field condition is satisfied and r far-field acoustic signals impinge on the array from distinct angles {θ 1 , θ 2 , . . . , θ r }, and M snapshots are collected by the array, the data matrix can thus be written as ., a 0 (θ r )] ∈ C N 0 ×r denotes the steering matrix of all the steering vectors, S ∈ C r×M is the signal matrix, N represents the additive Gaussian noise matrix, and the size of X 0 is N 0 × M. Moreover, we use M 0 = A 0 (θ)S ∈ C N 0 ×M to denote the measurement matrix, and thus the array data matrix can be expressed as X 0 = M 0 + N.
In array signal processing, the number of sources is usually much smaller than that of the array size, and we have r < min{N 0 , M}. Therefore, in the noiseless case, the measurement matrix M 0 is usually low-rank, and its rank corresponds to the number of sources [27,28]. Note that the low-rank characteristics of M 0 could be destroyed when M 0 is corrupted by the noise, which is very common in ocean environments. For illustration, in this work, we assume that the UUVS forms a large-aperture virtual linear array. More specifically, the array consists of U subarrays, and there is a certain distance between two adjacent UUVs. In other words, the UUVS array can be regarded as a sparse linear array. Unlike the traditional DOA estimation with sparse arrays, herein the array is considered to be a large-aperture array with virtual sensors (by assuming that the sensors between the UUVs fail to work properly), which cause unobserved entries in the data matrix, as demonstrated in Figure 3. As mentioned above, in general, the noiseless array data matrix has a low-rank property, which is of great importance for data recovery by using matrix completion (MC) [29][30][31][32]. Therefore, if we can recover the unobserved data (gray blocks in Figure 3), then a virtual array with more sensors than those of the actual array can be achieved. Accordingly, a better DOA estimation performance is expected.

Preliminaries of MC
The compressed sensing (CS) theory [33][34][35] pointed out that, under certain conditions, the signal can be under-sampled at a lower sampling rate than the Nyquist sampling frequency, and the original signal can be recovered through a small amount of sampled data. With rational utilization of sparsity, the performance of sonar systems can be greatly improved [36]. According to the CS theory, the sparsity of the targets in the spatial domain can be used to convert the received array data into a sparse representation problem, which can be adopted for spatial spectrum estimation [27]. The theory of MC was developed on the heels of CS [37], which indicates that a problem has nothing to do with a large number of entries but actually only depends on a small number of entries. If the position and number of observed entries meet certain conditions, the MC technique can correctly recover the unknown entries with a high probability [30][31][32].
From the MC theory, it is known that the entire matrix is reconstructed by finding the unknowns of the data matrix such that the rank of the matrix is minimized. This is different from those problems (e.g., [38]) using the sparsity property to recover a sparse variable. Mathematically, let X be a low-rank matrix of which only a sampled set of entries [X] ij , (i, j) ∈ Ω are known, then the MC problem can be written as [29] minimize rank(M) subject to P Ω (M) = P Ω (X) (2) where rank(·) denotes the matrix rank, and P Ω is a sampling operator with the (i, t)-th ∈ Ω, and 0 otherwise. An illustration for this is given in Figure 4, for which the diagonal entries are unobserved after sampling, but can be correctly recovered by using MC. The MC technology has been used in the recommendation system [39] to record the users' scoring data of product and make recommendations to them by analyzing their historical record. Moreover, it has been used in many engineering fields, such as image processing [40], system recognition [41], sensor networks [42], spectrum sensing [43], sparse channel estimation [44], multimedia coding communication [45], and genotype imputation [46]. In addition, since the signal generally satisfies the sparsity or low rank in the spatial and frequency domain, MC has also been successfully applied in array signal processing [47]. With accurate recovery of the data matrix and reducing hardware requirements of the array system, MC has shown promising usage in array failure [28], noise reduction [48], and down sampling [49]. Furthermore, the variant of MC, i.e., tensor completion, has also been exploited for large-scale array systems [50].

UUVS Array Data Pre-Processing
It can be seen from Figure 3 that the data matrix of a UUVS large-aperture uniform linear virtual array with N sensors has the following form: where x 1 , x 2 , · · · are the columns of the actual data matrix X T 0 , × denotes the entries of virtual sensors (which are regarded as unobserved entries), and (·) T is the transpose operator. However, when the entries of a certain row (column) of a matrix are completely unobserved, those entries cannot be directly recovered by solving the problem (2). An effective way to deal with this issue is reshaping the matrix before data recovery. In this work, we employ the Hankel matrix transformation based MC [28] to reshape the matrix and recover unobserved entries.
For Hankel matrix transformation, let us define a linear operator as H n 1 : X ∈ C M×N → X h ∈ C Mn 1 ×n 2 , then the original matrix X is mapped into the corresponding Hankel matrix as where x i is the ith column of X, and n 1 + n 2 = N + 1. It can be found that, even if some columns of the original data matrix are unobserved, those entries are rearranged in the reconstructed Hankel matrix, such that the case of unobserved columns will no longer appear. For a uniform linear array, we can choose n 1 = N/2 , and the data matrix after Hankel transformation is still low-rank. This allows us to reconstruct the UUVS virtual array data for DOA estimation. Applying the MC to X h yields a recovered matrix M h , which should be reshaped back to the same size as M by using the inverse Hankel transformation H † : where w j is the number of elements in the j-th anti-diagonal of an n 1 × n 2 matrix [28].
To better understand the transformation and data recovery, an example is given below. Suppose that the original data matrix M is given by Let Ω be the set of indexes of observed entries, then the observed matrix can be expressed as The corresponding Hankel structure matrix is where indexes of unobserved entries become {(1, 2), (2, 2), (3, 2), (4, 1), (5, 1), (6, 1)}, and the indexes of remaining entries form the set including the observed entries. The new set is denoted as Ω. More exactly, we have X h = P Ω (M h ) = P Ω (H(M)). By applying the MC, the unobserved entries in (8) can be recovered, and the matrix is reshaped according to (5).

UUVS Array Data Recovery via Structured MC
According to the aforementioned analysis of Hankel matrix transformation, in order to recover the low-rank matrix from the array data matrix, we formulate the structured MC problem as minimize rank(H(M)) subject to H(P Ω (X)) − H(P Ω (M)) 2 2 ≤ where · 2 denotes 2 norm and a tolerance parameter is introduced to take the noise in the array data matrix into account. Particularly, in the noiseless case, we can set = 0. According to the Hankel transformation described above, the problem (9) can be written as which is a rank minimization problem known to be Np-hard [30][31][32]. To solve the rank minimization problem (10), we relax the rank function to the kernel norm of matrix and the optimization problem is recast as where · * denotes the nuclear norm. It should be noted that the above problem formulation only takes the unobserved entries into account. In the case that the certain observed entries are corrupted with outliers, the solution may be unsatisfactory. To deal with this issue, a sparse matrix S is introduced to the data matrix M as M + S. Moreover, the sampled Hankel transformed matrix satisfies Z h = P Ω (Z h ) and it is also sparse. Therefore, the problem of recovering the low-rank matrix with unobserved and corrupted entries can be formulated as where · F and · 0 denote Frobenius norm and 0 norm, respectively, and s represents the maximal number of non-zero entries in Z h .
In order to solve the non-convex optimization problem (12), the concept of robust matrix completion (RMC) [51] approach is employed in this work. Briefly speaking, the RMC contains two iterative loops: in the k-th iteration of the outer loop, X h is decomposed into a matrix with rank k and a sparse matrix, which will be served as the initialization of the (k + 1)-th iteration. In the t-th iteration of the inner loop, the sparse matrix Z t+1 h and the matrix M t+1 h both with rank k are respectively updated based on the Z t h and M t h . Next, M t+1 h is computed with a gradient descent step, and then projected back to the set of rank-k matrices, i.e., where | Ω| denotes the size of set Ω, P k (A) stands for the projection of A onto the set of rank-k matrices, and the singular-value decomposition (SVD) can be adopted for this purpose. Specifically, let the SVD of A as A = ∑ i=1 σ i u i v T i with singular values σ i are arranged in descending order, then Next, Z t+1 h is determined by using the residual difference between X h and M t+1 h with a hard-thresholding operator as

DOA Estimation Based on Recovered Data Matrix
Having the recovered matrix M h by using the above RMC method, the estimate of the noiseless UUVS data matrix of the large-aperture virtual array can be obtaining by the inverse Hankel transformation (5) as On this basis, the traditional high-resolution DOA estimation approaches can be applied to determine the angles of acoustic sources. As a representative approach, the classic MUSIC algorithm is employed in this work. Nevertheless, other methods such as ESPRIT can also applied. Let the eigenvalue decomposition of M be where U s and U n denote the signal subspace and noise subspace, respectively, Σ s and Σ n contains the corresponding eigenvalues. The spatial spectrum for DOA estimation can thus be given by where a(θ) is the steering vector of the UUVS large-aperture virtual array.

Numerical Results
As shown in Figure 5, the scenario of a single line abreast composed of three UUVs is considered, and each UUV carries a 3-element hydrophone line array. Meanwhile, we assume that a fixed formation spacing is maintained among adjacent UUVs. Consequently, the UUVS with a single line abreast has multiple short linear arrays. As the size of a single UUV is limited, the aperture of hydrophone array carried by itself is small, which results in the low spatial gain and resolution. In this section, all the hydrophones in the single line abreast is regarded as a linear sparse array with N 0 = 9 physical sensors. With the MC technique, we can extend the array to a large-aperture array with N = 15 sensors (the six sensors between the UUVs are virtual.) More specifically, the virtual array sensors' positions are {4, 5, 6, 10, 11, 12}, and the UUV array elements and virtual array elements form a 15-element virtual linear array with inter-element spacing being half-wavelength. Hence, for the virtual array, the failure percentage of the array elements is 40%. The signals are assumed to be monochromatic waves with frequency 750 Hz and the signal-to-noise ratio (SNR) defined as SNR = 10 log 10 σ 2 s /σ 2 n dB is set to be 20 dB. Here, σ 2 s and σ 2 n are the variances of signal and noise, respectively. Additive Gaussian noise is considered. We take M = 500 snapshots and conduct Q = 200 Monte Carlo experiments. The performance of the proposed DOA estimation method via UUVS is examined under various parameter settings. The root-mean-square-error (RMSE) of DOA estimation is where θ q,p denotes the DOA estimate of the q-th signal in the p-th experiment. Meanwhile, different MC algorithms, including singular value thresholding (SVT) [52], OPTSPACE [53] and 2 approximation [54] are used to compare the performance of acoustic target DOA estimation under the UUVS scheme. The performance of data recovery algorithm is measured in terms of the normalized recovery error (NRE), which is defined as NRE( M) = E{ M − M 2 F / M 2 F }.

Single Target
In this example, the number of target signals is r = 1, and the incident angle is 50 • . Figure 6 shows the data NREs achieved by different MC algorithms. It is seen that the NRE performance in the ascending order is: SVT, OPTSPACE, 2 approximation, and RMC.  Figure 7 compares the MUSIC spectra for DOA estimation. It can be seen that the array of neither a single UUV nor multiple UUVs can satisfactorily estimate the target angle. All four of the methods with structured MC can work well and successfully estimate the incident angle of target. However, their peak sharpnesses are somewhat different. Specifically, the 2 approximation and RMC algorithms have the best DOA estimation performance. It is also worth noting that the RMSE(θ) of DOA estimation with a single UUV is 4.3513 • , and that with multiple UUVs is 0.03125 • . When the received data are recovered via the structured MC, the RMSE(θ) of SVT, OPTSPACE, 2 approximation and RMC are all zero.

Two Targets
Now, we assume r = 2, θ 1 = 25 • and θ 2 = 50 • . Figure 8 shows the data NRE achieved by different structured MC algorithms. In this case, the NRE from high to low is: SVT, OPTSPACE, 2 approximation, and RMC. Figure 9 shows the DOA estimation spectra of two targets based on the proposed methods. It is seen that the RMC algorithm exhibits the best DOA estimation performance among all the structured MC methods. The RMSE(θ) of DOA estimation with a single UUV is 4.0456 • , and that with multiple UUVs is 4.5781 • . For the received data reconstructed via the structured MC, the RMSE(θ) of OPTSPACE is 6.25e−04 • , and those of SVT, 2 approximation and RMC are all zero. Compared with a single target, the data recovery error and DOA estimation results of the two targets are different, which implies that the rank of matrix has a significant influence on data recovery and DOA estimation.

Large Number of Targets
In this example, the number of target signals is set to r = 5, and the incident angles are {10 • , 25 • , 40 • , 50 • , 70 • }, respectively. Figure 10 investigates the data NRE of different structured MC algorithms. Figure 11 plots the DOA estimation spectra of five targets. We can observe that the array of neither a single UUV nor multiple UUVs can successfully estimate the targets. At the same time, the methods of structured MC with SVT, OPTSPACE, and 2 approximation are unable to estimate the targets, but the RMC can still successfully estimate the incident angles of five targets. However, when the number of signal targets is 6, none of the aforementioned methods can successfully estimate the incident angles. The advantages of UUVS can be enhanced based on the proposed methods. It is also worth noting that, compared with a single UUV, the array aperture has been significantly increased, and compared with multiple UUVs, the array gain of the equivalent virtual full array can be increased through data recovery. Therefore, the estimatable target number has been significantly increased, compared with the traditional methods. Additionally, the RMSE(θ) of DOA estimation with multiple UUVs is 0.2425 • , and the RMSE(θ) of SVT, OPTSPACE, 2 approximation and RMC are 0.00125 • , 0.0165 • , 0.00325 • and 0, respectively. In conclusion, the performance of RMC is suitable for underwater target DOA estimation via UUVS.

Angular Resolution
Following [55,56], the discriminant threshold formula of angular resolution is defined as follows: Υ(θ 1 , θ 2 ) P (θ m ) − 1/2{P (θ 1 ) + P (θ 2 )}, (19) where P (θ) = 1/P(θ), and P(θ) denotes the MUSIC spectrum, and P (θ) is the corresponding null spectrum. In addition, θ 1 and θ 2 are the incident angles of two targets, and θ m = (θ 1 + θ 2 )/2 defines the central angle. The angular interval of targets is represented as ∆θ, so we have θ 1 = θ m − ∆θ/2, and θ 2 = θ m + ∆θ/2. Additionally, we use P res {Υ > 0} to denote the successful probability of DOA resolution, i.e., those two incident targets are successfully distinguished; otherwise, they cannot be successfully distinguished. When θ m is 45 • , θ 1 and θ 2 are, respectively, 44.5 • and 45.5 • , so the angular interval is 1 • . Figure 12 shows the data NRE achieved by different structured MC algorithms when the target angular interval is 1 • . The performance in the ascending order is: SVT, OPTSPACE, 2 approximation, and RMC. Furthermore, it is observed that the RMC can successfully distinguish two incident targets when the angular interval is 1 • , as shown in Figures 13 and 14. The success probability of DOA resolution for different angular intervals is shown in Figure 15. We can observe that, when the angular intervals are 1 • and 2 • , the SVT cannot successfully distinguish the adjacent targets; when the angular interval is larger than 2 • , the SVT can successfully distinguish adjacent targets. Moreover, the OPTSPACE cannot distinguish targets with close angular intervals, and it has the worst performance. Note that, as the angular interval increases, the probability of 2 approximation that can successfully distinguish the adjacent targets becomes higher, and the resolution performance gradually improves. Finally, RMC exhibits the best resolution performance, as it can successfully distinguish adjacent targets with different angular intervals.

Conclusions
UUVS is an important development trend of marine equipment and has outstanding potential for underwater acoustic target localization. This paper presents a DOA method of underwater acoustic sources with UUVS via structured MC. The mathematical model of UUVS array composed of multiple subarrays is described and an effective algorithm of structured MC is introduced. Specifically, a virtual full array is modeled based on multiple UUVs, and the data of virtual array elements are recovery by MC based on Hankel transformation. As the noiseless low-rank matrix of the equivalent large-aperture virtual array is recovered, the DOA estimation for underwater target can be obtained by the MUSIC algorithm. Simulation results demonstrate the performance of the proposed method. Compared with the traditional single array method, the proposed method with UUVS has various advantages including higher accuracy of DOA estimation, higher resolution, and more estimatable targets.