Partial Angular Sparse Representation Based DOA Estimation Using Sparse Separate Nested Acoustic Vector Sensor Array

In this paper, the issue of direction of arrival (DOA) estimation is discussed, and a partial angular sparse representation (SR)-based method using a sparse separate nested acoustic vector sensor (SSN-AVS) array is developed. Traditional AVS array is improved by separating the pressure sensor array and velocity sensor array into two different sparse array geometries with nested relationship. This improved array geometry can achieve large degrees of freedom (DOF) after the extended vectorization of the cross-covariance matrix, and only partial SR of the angle is required by exploiting the cyclic phase ambiguity caused by the large inter-element spacing of the virtual array. Joint sparse recovery is developed to amend the grid offset and unitary transformation is utilized to transform the complex atoms into real-valued ones. After sparse recovery, the sparse vector can simultaneously provide high-resolution but ambiguous angle estimation and unambiguous reference angle estimation embedded in the AVS array, and they are combined to obtain unique and high-resolution DOA estimation. Compared to other state-of-the-art DOA estimation methods using the AVS array, the proposed algorithm can provide better DOA estimation performance while requiring lower complexity. Multiple simulation results verify the effectiveness of the approach.


Introduction
Sensor arrays can utilize signals from multiple paths to overcome fading effect and enhance system capacity, so they have found wide application in many fields [1][2][3][4]. As a key issue for sensor array signal processing, direction of arrival (DOA) estimation has attracted lots of attention in many systems, e.g., radar [5,6], sonar [7], and wireless communication [8,9], to name a few. Compared to the conventional scalar sensor array, which only measures signal pressure information, the acoustic vector sensor (AVS) array provides a vector output, which contains both pressure and velocity information, so it can further extend array aperture and strengthen parameter identifiability [10]. Furthermore, the AVS array itself can provide unambiguous reference DOA estimation, which enables the inter-element spacing of the array to be larger than half-wavelength to enhance spatial resolution [11]. Since its introduction, DOA estimation using AVS array has attracted growing interest [11,12]. Many effective DOA estimation methods using AVS array have been developed, such as the estimation of signal parameters via rotational invariance technique (ESPRIT)-based method [11], Capon-based method [12], sensors are extracted and arranged along the negative side with inter-element spacing being MNd. The pressure sensor and velocity sensors are no longer co-located but separated into two different geometries with nested relationship. Assuming that there are K far-field uncorrelated source signals impinging on the proposed array, then the output of the pressure sensor array is formulated as Similarly, the outputs of the velocity vector sensor array are expressed as   (5) According to Equation (2) and Equation (5), the steering vectors of the pressure sensor array and velocity sensor array satisfy nested relationship, which will greatly increase the DOF in the co-array domain shown in the next section. It is also noted that the proposed SSN-AVS array just separates the pressure array and velocity array, it does not add extra AVS, which can be observed from Figure 1, so the total AVS number is still M. Assuming that there are K far-field uncorrelated source signals impinging on the proposed array, then the output of the pressure sensor array is formulated as

Virtual Array in the Co-Array Domain
where s(t) = [s 1 (t), s 2 (t), · · · , s K (t)] T is a vector containing all the K source signals, n p (t) is the additive white Gaussian noise vector, and A p = [a p (θ 1 ), a p (θ 2 ), · · · , a p (θ K )] denotes the direction matrix, which depends on the geometry of pressure sensor array. The steering vectors a p (θ k ), k = 1, . . . , K are expressed as Similarly, the outputs of the velocity vector sensor array are expressed as where Φ 1 = diag(sin θ 1 , · · · , sin θ K ) and Φ 2 = diag(cos θ 1 , · · · , cos θ K ) are two diagonal matrices reflecting different components along different axes. n v1 (t) and n v2 (t) are additive white Gaussian noise vectors.
is the direction matrix of the velocity sensor array, and the steering vectors are expressed as According to Equation (2) and Equation (5), the steering vectors of the pressure sensor array and velocity sensor array satisfy nested relationship, which will greatly increase the DOF in the co-array domain shown in the next section. It is also noted that the proposed SSN-AVS array just separates the pressure array and velocity array, it does not add extra AVS, which can be observed from Figure 1, so the total AVS number is still M.
Thereafter, to obtain the virtual array in the co-array domain, we vectorize the cross-covariance matrix as T is a vector, and can be regarded as a new signal vector in the co-array domain.
can be regarded as a new direction matrix corresponding to a M 2 -element virtual array, whose steering vectors are expressed as Based on Equation (2) and Equation (5), now an M 2 -element virtual array with inter-element spacing Nd can be achieved. Compared to the original physical array in Figure 1, the virtual array shows much larger DOF. Additionally, the new signal vector p 1 is a real-valued one, which means that the array aperture can be furtherly extended by combining the conjugate data of the vector in Equation (7). Construct the extended vector as where y vp1 (1 : M 2 − 1) means to extract the first (M 2 −1) elements of y vp1 , and A pv (1 : M 2 − 1, ) means the first (M 2 − 1) rows of A pv . Here, we extract the first (M 2 − 1) elements of y vp1 to eliminate the overlapped origin element in the extended virtual array.
T is the extended direction matrix, whose steering vectors a(θ k ), k = 1, . . . , K are expressed as Based on Equation (9), now a (2M 2 − 1)-element virtual array with steering vectors being a(θ k ) can be obtained, and the large DOF can be expected to improve the parameter estimation performance and system capacity. If the traditional array geometry in Figure 1a is used, then the virtual steering vector in Equation (8) will have many redundant elements due to the same geometry being used by the pressure array and velocity array. Additionally, the aperture extension in Equation (9) will be unavailable, and the final achievable DOF will be (2M − 1), which is much smaller than what the proposed geometry achieved.
It is shown that the steering vector in Equation (10) satisfies conjugate symmetrical property, so unitary transformation can be employed to transform the complex vector into real-valued one [30], which can effectively reduce the computation complexity. The unitary matrix is defined as Then the extended steering vector can be transformed into a real-valued vector According to Equation (9) and Equation (12), the real-valued data is obtained via Based on x v1 (t), x p (t) and the steps from Equation (6) to Equation (13), now the real-valued data y rz1 is obtained. Similarly, based on x v2 (t) and x p (t), whose cross covariance matrix is where Combining y rz1 and y rz2 yields where p = [e jθ 1 σ 2 1 , . . . , e jθ K σ 2 K ] contains unambiguous DOA information, which is embedded in the AVS array (generated from Φ 1 and Φ 2 ). Now a virtual array with output being y z is obtained, and the direction matrix A r is real-valued, and sparse representation-based DOA estimation can be developed.

Partial Angular Sparse Representation
Suppose that the angle grid of interest is θ 1 , θ 2 , · · · , θ P (P K), which contains the true DOAs. Then construct the dictionary Ω = [a r ( θ 1 ), . . . , a r ( θ P )], which must contain the columns of A r , so Equation (15) can be written as where ρ is a P-element vector, which only has K non-zero elements from p, and the positions of the non-zero elements give the DOA estimations. However, according to Equation (10), the inter-element spacing of the (2M 2 − 1)-element virtual array is Nd, which means phase ambiguities in this array. Let z p = e jNπ sin θ p denotes the phase difference between adjacent elements corresponding to angle θ p , then there are another (N − 1) angles θ p,n , n = 2, . . . , N satisfying e jNπ sin θ p,n = z p , n = 2, . . . , N The phase ambiguity shown in Equation (17) will cause multiple false angle estimations. However, we can in turn exploit the ambiguity to shrink the size of the dictionary and focus on one specific cycle. It is derived from Equation (17) that the N angles satisfy where n v is an integer ensuring sin θ p,n located in the range [−1, 1]. According to Equation (18), the N solutions (sin θ p and sin θ p,n , n = 2, . . . , N) are uniformly distributed within the range [−1, 1] with adjacent distance being 2/N. Figure 2 shows an example when sin θ p = 0.5 and N = 3, and other two ambiguous solutions are sin θ p,2 = −1/6 and sin θ p,3 = −5/6, respectively. If we divide the range [−1, 1] into N cycles with each width being 2/N, then one cycle provides only one solution. Therefore, we can choose one specific cycle as a representative to construct the dictionary, and recover other solutions based on their uniform distribution after sparse vector recovery. Without loss of generality, we choose range [−1/N, 1/N], whose corresponding angular range and dictionary are [arcsin(−1/N), arcsin(1/N)] and Ω rep , respectively. Then the sparse representation can be expressed as where the elements in ρ ∈ C P×1 corresponding to the representative solutions are the same as those in p, and the others are zero. For example, as Figure 2 shows, the true solution 0.5 is represented by −1/6 in Ω rep .  (18) where v n is an integer ensuring where the elements in Now the sparse representation based on partial angular domain is conducted, and the phase ambiguity can be avoided. However, Equation (19) is established based on the assumption that the representative angles are right located in the pre-defined grid, but the angles are very likely to lie off the grid, no matter how fine the grid is defined. The grid mismatch will degrade the sparse recovery performance. Consequently, an improved sparse presentation framework should be developed to tackle the off-grid problem.

Joint Sparse Representation Framework
In this section, we take the off-grid angles into account, and utilize the joint sparse representation framework to enhance the robustness to grid mismatch.
Denote the new grid within the partial angular range Then, Equation (19) has a robust form Now the sparse representation based on partial angular domain is conducted, and the phase ambiguity can be avoided. However, Equation (19) is established based on the assumption that the representative angles are right located in the pre-defined grid, but the angles are very likely to lie off the grid, no matter how fine the grid is defined. The grid mismatch will degrade the sparse recovery performance. Consequently, an improved sparse presentation framework should be developed to tackle the off-grid problem.

Joint Sparse Representation Framework
In this section, we take the off-grid angles into account, and utilize the joint sparse representation framework to enhance the robustness to grid mismatch.
Denote the new grid within the partial angular range [arcsin(−1/N), arcsin(1/N)] as θ 1 , θ 2 , · · · , θ Q , which is uniformly sampled. Then the true representative angle θ k can be represented by a nearest grid angle θ q,k plus an offset α k (maybe zero when the true angle is right on the gird). Based on first-order Taylor expansion around the predefined grid, the true steering vector can be approximately expressed as [28,31] a(θ k ) ≈ a( θ q,k ) + ∂a( θ q,k ) Then, Equation (19) has a robust form where Ω sub = ∂a( θ 1 ) denotes the derivatives on the sampling grid, Λ = diag(β) is a diagonal matrix containing the offsets, and So β is also a sparse vector with K non-zero elements. Let ω = Λρ, then ω and ρ are joint sparse. Equation (21) is rewritten as Equation (23) is a joint sparse recovery problem, and joint orthogonal matching pursuit (JOMP) [31] can be utilized to recover the sparse vectors ρ and ω. The basic steps of JOMP are shown in Table 1. It is noted that the grid can be coarse due to offset compensation afterwards [28], so the sparse recovery complexity can be reduced.

Ambiguity Elimination and DOA Estimation
After sparse recovery using JOMP, the estimations of ρ and ω can be obtained. The positions of non-zero elements in ρ give the grid angles θ q,k , k = 1, . . . , K, which are nearest to the representative angles. Meanwhile, the offsets are obtained via From Equation (24), the offsets α k , k = 1, . . . , K can be acquired, then the angles estimated in the representative range are The angles obtained in Equation (25) are representative angles. Finally, we need to resolve the ambiguity problem in the angle estimation. According to Equation (18), from one specific angle θ k , we can recover a total of N angles.
To determine the true angle in the N angles, we need an unambiguous reference angle estimation, which is actually embedded in the sparse vector ρ, whose non-zero elements give the estimation of p = [e jθ 1 σ 2 1 , . . . , e jθ K σ 2 K ]. Then the unambiguous reference angle is obtained via θ k = angle(p(k)), k = 1, . . . , K By finding the angle nearest to θ k in the N angles, the true DOA estimation can be determined. It should be noted that the ambiguous angle and unambiguous reference angle are automatically paired due to the sparse vector.

Remarks and Summary
Remark 1: In practice, the cross-covariance matrix in Equation (6) is estimated via finite snapshots where T denotes the number of snapshots. Remark 2: One-dimensional angles are assumed in this paper, which can find its applications in some situations where the signals are in the same plane with the arrays, such as sea-surface detection, satellite-to-satellite location, ground signal finding and so on.
Remark 3: For DOA estimation methods, including the proposed method, the source number K is generally assumed to be known a priori. The source number estimation is another important issue for array signal processing [32].
The major steps of the proposed algorithm are: 1.
Construct the cross-covariance matrix between the pressor sensor array and velocity sensor array via Equation (27); 2.
Vectorize the cross-covariance matrix to obtain a vector and combine its conjugate to form a (2M 2 − 1)-element virtual array with inter-element spacing being Nd; 3.
Employ unitary transformation to transform the complex direction matrix into a real-valued one and form the final output via Equation (15) After sparse recovery based on JOMP, obtain the grid and offset estimations and combine them to obtain the representative angle estimation using Equation (25); 6.
Recover all angle estimations based on the representative angle, and determine the final DOA estimation based on the unambiguous angle embedded in the sparse vector.
Based on the steps shown above, the proposed method requires neither EVD nor peak searching, and only requires a real-valued dictionary covering the partial angular domain, so the main complexity of the proposed method lies in cross-covariance matrix construction and joint sparse recovery. The theoretical complexity (measured by the number of complex multiplications) of the proposed method is shown in Table 2, where the main complexities of the Successive MUSIC [17], tensor-based method [25] and SR-based method [26] are also presented for comparison. The same AVS numbers are used, and n 1 denotes the local search steps used by Successive MUSIC, M 1 = M 2 /4 + M/2, n 2 denotes the global search steps used by tensor-based method, n 3 denotes the dictionary size of the SR-based method, n 4 denotes the dictionary size of the proposed method. A typical setting is: M = 4, T = 200, K = 2, n 1 = 60, n 2 = 90, n 3 = 90, n 4 = 60. It is indicated in Table 2 that the proposed method costs less computation resources than other methods.

Simulation Results
In the simulations, the proposed SSN-AVS array is configured as M = 4 and N = 2, which means there are totally 4 AVSs. For the other methods, the successive MUSIC adopts traditional array shown in Figure 1a, the tensor-based method adopts nested AVS array and the SR-based method uses separated nested AVS array. The total AVS numbers used are the same for fair comparison. Assume there are two uncorrelated signals with DOAs being θ 1 = 20.76 • and θ 2 = 45.34 • , respectively. T = 200 snapshots are collected, and the root mean square error (RMSE) of the angle estimation is defined below to measure the DOA estimation performance whereθ k,l denotes the estimation of θ k of the l-th Monte Carlo trial, whose total number is L = 500. Figure 3 shows the DOA estimation results of the proposed algorithm over 100 trials when SNR = 0 dB, and it is shown that both the two DOAs can be accurately estimated. In the simulations, the proposed SSN-AVS array is configured as M = 4 and N = 2, which means there are totally 4 AVSs. For the other methods, the successive MUSIC adopts traditional array shown in Figure 1a, the tensor-based method adopts nested AVS array and the SR-based method uses separated nested AVS array. The total AVS numbers used are the same for fair comparison. Assume there are two uncorrelated signals with DOAs being 1  snapshots are collected, and the root mean square error (RMSE) of the angle estimation is defined below to measure the DOA estimation performance where ,kl  denotes the estimation of k  of the l-th Monte Carlo trial, whose total number is L = 500. Figure 3 shows the DOA estimation results of the proposed algorithm over 100 trials when SNR = 0 dB, and it is shown that both the two DOAs can be accurately estimated. The DOA estimation performance comparisons between the successive MUSIC [17], tensorbased method [25], SR-based method [26] and the proposed algorithm versus SNR and snapshot number are shown in Figure 4 and Figure 5, respectively. The grid interval used by the successive MUSIC, tensor-based method and SR-based method is 0.1°, while the grid interval of the proposed method is 1°, which can save much complexity in the sparse recovery. It is shown in Figure 4 and Figure 5 that the proposed method achieves better DOA estimation performance than the other three methods. One reason is the large DOF and aperture generated by the proposed SSN-AVS array, and the other reason is that the proposed method can amend the grid offsets caused by the off-grid sources, which lead to the performance degradations of the other three methods, despite a finer grid they used. To verify the two reasons that lead to the performance improvement, we test the proposed method under other two configurations: 1. uses traditional array geometry, 2. uses the proposed geometry but utilizes general SR instead of joint SR in Section 3.3. The DOA estimation performance comparison between these three situations has been shown in Figure 6, which indicates that the both the array geometry and joint SR contribute to the performance improvement, especially the array geometry. As explained in Section 3.3, the aperture extension is unavailable and the final achievable DOF will be greatly reduced if the traditional array geometry is used. The DOA estimation performance comparisons between the successive MUSIC [17], tensor-based method [25], SR-based method [26] and the proposed algorithm versus SNR and snapshot number are shown in Figures 4 and 5, respectively. The grid interval used by the successive MUSIC, tensor-based method and SR-based method is 0.1 • , while the grid interval of the proposed method is 1 • , which can save much complexity in the sparse recovery. It is shown in Figures 4 and 5 that the proposed method achieves better DOA estimation performance than the other three methods. One reason is the large DOF and aperture generated by the proposed SSN-AVS array, and the other reason is that the proposed method can amend the grid offsets caused by the off-grid sources, which lead to the performance degradations of the other three methods, despite a finer grid they used. To verify the two reasons that lead to the performance improvement, we test the proposed method under other two configurations: 1. uses traditional array geometry, 2. uses the proposed geometry but utilizes general SR instead of joint SR in Section 3.3. The DOA estimation performance comparison between these three situations has been shown in Figure 6, which indicates that the both the array geometry and joint SR contribute to the performance improvement, especially the array geometry. As explained in Section 3.3, the aperture extension is unavailable and the final achievable DOF will be greatly reduced if the traditional array geometry is used. Sensors 2018, 17, x FOR PEER REVIEW 10 of 15            In Figures 7 and 8, two closely spaced sources with DOAs respectively being θ 1 = 20.76 • and θ 2 = 24.34 • .34 • are adopted to test the resolutions of the algorithms. When SNR = 5 dB, subspace-based methods including the tensor-based method and successive MUSIC fail to identify the two closely spaced sources, while the SR-based method and the proposed method still work well due to the super-resolution property of sparse recovery technique. When SNR gets lower in Figure 8, the SR-based method fails, and the proposed method can approximately identify the two sources (with visible deviations). Consequently, the proposed method shows finer angular resolution than the other methods. To further evaluate the resolution performances of the methods, Figure 9 shows the DOA estimation performance comparison versus angular separation when SNR = 5 dB. The reference angle is set as θ 1 = 40.36 • , and the second angle is θ 2 = θ 2 + ∆θ, where ∆θ denotes the angular separation varying among the range [2 • , 6 • ]. It is shown that the subspace-based methods cannot achieve effective results within the given range, and the SR-based method and the proposed method work better, especially when ∆θ > 4 • (can clearly identify the two sources). Meanwhile, the proposed method always achieves the best DOA estimation performance versus the variation of the angular separation. are adopted to test the resolutions of the algorithms. When SNR = 5 dB, subspacebased methods including the tensor-based method and successive MUSIC fail to identify the two closely spaced sources, while the SR-based method and the proposed method still work well due to the super-resolution property of sparse recovery technique. When SNR gets lower in Figure 8, the SR-based method fails, and the proposed method can approximately identify the two sources (with visible deviations). Consequently, the proposed method shows finer angular resolution than the other methods. To further evaluate the resolution performances of the methods, Figure 9 shows the DOA estimation performance comparison versus angular separation when SNR = 5 dB. The reference angle is set as 1 40.36

 =
, and the second angle is 22    = +  , where   denotes the angular separation varying among the range [2 , 6 ] . It is shown that the subspace-based methods cannot achieve effective results within the given range, and the SR-based method and the proposed method work better, especially when 4   (can clearly identify the two sources). Meanwhile, the proposed method always achieves the best DOA estimation performance versus the variation of the angular separation.   Figure 10 shows the DOA estimation performance comparison with non-uniform spatial noise. Compared to Figure 4, both the tensor-based method and successive MUSIC method have performance degradations. The subspace-based methods are sensitive to nonuniform noise due to the permeation between the signal and noise subspaces. Both the SR-based method and the proposed method is robust to nonuniform noise due to the usage of cross covariance matrix, which can effectively restrain the power of nonuniform noise.    Figure 10 shows the DOA estimation performance comparison with non-uniform spatial noise. Compared to Figure 4, both the tensor-based method and successive MUSIC method have performance degradations. The subspace-based methods are sensitive to nonuniform noise due to the permeation between the signal and noise subspaces. Both the SR-based method and the proposed method is robust to nonuniform noise due to the usage of cross covariance matrix, which can effectively restrain the power of nonuniform noise.   Figure 10 shows the DOA estimation performance comparison with non-uniform spatial noise. Compared to Figure 4, both the tensor-based method and successive MUSIC method have performance degradations. The subspace-based methods are sensitive to nonuniform noise due to the permeation between the signal and noise subspaces. Both the SR-based method and the proposed method is robust to nonuniform noise due to the usage of cross covariance matrix, which can effectively restrain the power of nonuniform noise. Sensors 2018, 17, x FOR PEER REVIEW 13 of 15 Figure 10. DOA estimation performance comparison with non-uniform spatial noise.

Conclusions
In this paper, the SSN-AVS array is designed and a corresponding partial angular SR-based DOA estimation method is proposed. Multiple analyses and simulations verify that the proposed method has the following advantages: