Particle-Velocity Coarray Augmentation for Direction Finding with Acoustic Vector Sensors

In this paper, the problem of passive direction finding is addressed using an acoustic vector sensor array (AVS), which may be deployed either in free space or near a reflecting boundary. Building upon the 4×1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4 \times 1$$\end{document} vector field measured by an AVS, the particle-velocity coarray augmentation (PVCA) is proposed to admit the underdetermined direction finding using the spatial difference coarray derived from the vectorization of the array covariance matrix. Unlike the widely used spatial coarray Toeplitz recovery technique, the PVCA is applicable to arbitrary array geometries and imposes no reduction of the spatial difference coarray aperture. For the array located at or near a reflecting boundary, the PVCA allows resolving up to 13 sources, while for the array located in free space, the PVCA can identify 9 sources at most. By applying to the systematically designed nonuniform arrays, such as coprime arrays and nested arrays, the PVCA can be coupled with the spatial smoothing technique to get the number of resolvable sources multiplied. Finally, the effectiveness of the PVCA is verified by numerical simulations.


Introduction
An acoustic vector sensor (AVS) is composed of one isotropic pressure sensor and up to three mutually orthogonal particle velocity sensors, all collocated at a single point in space [23]. The pressure sensor measures the acoustic pressure field as a scalar, while each velocity sensor measures the acoustic particle velocity in one of the corresponding Cartesian axes. An AVS can thus provide a 4 × 1 vector measurement of the acoustic field at an arbitrary spatial location [38]. In the context of passive array processing, an AVS array would offer improved performance over the pressure sensor array of the same aperture [9]. Acoustic vector sensors are commercially available from Microflown (http://www.microflown-avisa.com/). They are versatile in handling a wide variety of applications such as detection, localization, and tracking in air acoustic [1,21], battlefield surveillance and classification [2,10], and underwater communications [31].
Within the forgoing mentioned applications, direction finding is an essential problem as it plays a pivotal role in the overall success of a specific task, e.g., target localization and tracking. In fact, direction finding with AVS arrays has already been extensively investigated in the past decades. Many advanced algorithms have been reported in the open literature [5, 12, 14-16, 19, 25, 30, 35, 36, 41]. These algorithms are generally designed with the use of specific but regular array geometries, e.g., linear [5,12,15,19], circular [30,41], rectangular [14,35,36], L-shaped [25], and cross-shaped [16]. They are, of course, unable to work properly for an arbitrary array configuration. The recently developed tensor decomposition algorithm [34] is applicable to arbitrary array configuration. However, this algorithm, as well as the algorithms as given in [5, 12, 14-16, 19, 25, 30, 35, 36, 41], requires an overdetermined presumption, i.e., the number of source signals is less than that of the sensors, and consequently, being inapplicable to handle underdetermined direction finding.
On the other hand, the algorithms mentioned in the last paragraph consider the scenario of free space propagation, which implicitly assumes that the array is deployed far away from any reflecting boundary. However, in applications like hull-mounted sonar and battlefield surveillance, sensor arrays may be deployed at or near reflecting boundaries. For these applications, the algorithms designed under the free space assumption would fail. By taking the reflecting boundary into account, the general signal measurement model is introduced in [11], where the Cramér-Rao bound (CRB) for direction finding is derived as well. Moreover, the vector field smoothing algorithm and the ESPRIT-like algorithm have been developed for azimuth-elevation direction finding. Details may refer to [32] and [37]. However, these two algorithms are also based on the overdetermined assumption and are inapplicable to the underdetermined scenarios, where the number of source signals exceeds that of the sensors.
To overcome the foregoing drawbacks, this paper is aimed at proposing a new scheme of what is called particle-velocity coarray augmentation (PVCA) for underdetermined direction finding. The key insight of the PVCA is that it provides a full-rank augmented matrix to exploit the entire spatial difference coarray aperture, which is derived from the vectorization of the array covariance matrix. Compared with the many existing techniques, the PVCA offers the following four new properties: -First, the PVCA is applicable to arbitrary array geometry and does not reduce the spatial coarray aperture. This flexibility is not offered by the commonly used coarray processing techniques, e.g., the spatial coarray Toeplitz recovery technique [20], which decreases the spatial coarray aperture and can only apply to uniform coarrays -Second, the PVCA is applicable to scenarios, where the sensor arrays are deployed either in free space or near a reflecting boundary. In contrast, this adaptability is not offered by the existing algorithms. In fact, some of these algorithms work only for free space applications [5, 12, 14-16, 19, 25, 30, 35, 36, 41], while the others work only for the scenarios including the reflecting boundaries [11,32,37] -Third, it will be proved that for the array located at or near a reflecting boundary, the PVCA allows to resolve up to 13 sources, while for the array located in free space, the PVCA can identify 9 sources at most. Therefore, the PVCA is able to deal with the underdetermined cases for identifying more sources than sensors -Fourth, the PVCA can be applied to the systematically designed nonuniform arrays, such as minimum redundancy arrays [22,26,40], coprime arrays [27,33,39], and nested arrays [8,13,24]. In this way, it can be coupled with the spatial smoothing technique to get the number of resolvable sources multiplied The notation used throughout the paper is given in Table 1.

Signal Model and Problem Formulation
To set up the measurement signal model, the steering vector of the spatially collocated AVS is first introduced. For a source signal emitting from direction (θ, φ), the 4 × 1 steering vector of the AVS, which is located at or near a reflecting boundary, is given by [11] c(θ, φ) where 0 < θ < π denotes the elevation angle measured from the positive z-axis, 0 < φ < 2π denotes the azimuth angle measured from the positive x-axis, R is the complex reflection coefficient, representing the amplitude and phase variations of the reflected signal, ϑ = 4π d sin θ , d is the vertical distance between the sensor and the reflecting boundary, and u sin θ cos φ, v sin θ sin φ, and w cos θ denote the direction-cosines along the x-axis, y-axis, and z-axis, respectively. In (1), the first entry denotes the response of the pressure sensor, while the remaining three entries denote the responses of the three velocity sensors, where each aligns with one of the Cartesian axes. When the AVS is deployed far away from the reflecting boundary, Eq.
For an L-element array, the L ×1 spatial steering vector that relates the propagation delays between sensors to the source direction (θ, φ) is denoted as where λ denotes the system wavelength, r [x , y , z ] T denotes the location of the th sensor, and p [u, v, w] T denotes the source directional vector, i.e., the unit length vector pointing from the Cartesian origin toward the source. From (1) and (2), for an L-element AVS array, the 4L × 1 steering vector with respect to a source signal at direction (θ, φ) is expressed as With a total of K source signals, the data sample received by the entire AVS array at time t can be written as the following 4L × 1 vector: where θ k and φ k denote, respectively, the elevation angle and azimuth angle of the kth source signal, s k (t) denotes the waveform of the kth signal, Q [q(θ 1 , φ 1 ), · · · , q(θ K , φ K )] denotes the L × K manifold matrix that is associated with the array geometry, [a(θ 1 , φ 1 ), · · · , a(θ K , φ K )] denotes the 4L × K AVS array manifold, s(t) [s 1 (t), · · · , s K (t)] T denotes the K × 1 signal vector, and n(t) denotes the additive noise vector of size 4L × 1.
The problem at hand is to estimate the directions {(θ k , φ k ), k = 1, · · · , K } from the collected samples {x(t)} T t=1 . To that end, the following assumptions are made: 1. the source directions are distinct from one another, i.e., (θ 1 , φ 1 ) = (θ 2 , φ 2 ) = · · · = (θ K , φ K ); 2. the number of sources K is known or was successfully detected; 3. the source signals are modeled as zero-mean spatial-temporal uncorrelated stochastic processes. Thus, the source covariance matrix is of a diagonal structure and denoted as R s E{s(t)s H (t)} = diag[σ 2 1 , · · · , σ 2 K ], where σ 2 k is the power of the kth source signal; 4. the noise is spatially and temporally independent, identically distributed zero-mean white complex Gaussian, with a covariance matrix R n E{n(t)n H (t)} = σ 2 n I 4L , where σ 2 n denotes the deterministic but unknown variance of the noise; 5. The source signals are statistically independent of the noise.

Derivation of the PVCA
Under the assumptions made in Sect. 2, the 4L ×4L covariance matrix of x(t) is given by By performing the vectorization operation on R, a vector of size 16L 2 × 1 can be formed and is expressed as which is often called coarray output. In (6), ( A * A) = ( Q C) * ( Q C) denotes the 16L 2 × K spatial-particle-velocity coarray manifold, β [σ 2 1 , · · · , σ 2 K ] T , and i vec(I 4L ). Before proceeding to the PVCA, the rows of the spatial-particlevelocity coarray manifold need to be rearranged to the form of (C * C) ( Q * Q). Mathematically, this can be accomplished by left multiplying a 16L 2 × 16L 2 matrix P to ( A * A), i.e., The compact formation of P is derived in Appendix A. In (7),C C * C denotes the 16 × K particle-velocity coarray manifold matrix andQ Q * Q denotes the L 2 × K spatial coarray manifold matrix. Note that the manifold matrixC is irrelevant to the array geometry, while the manifold matrixQ is dependent only on the array geometry. This property forms the "core foundation" of the PVCA. Now, left multiplying r with P leads to a row-exchanged version of r, which can be expressed asr where e m is a vector of size m 2 × 1 with its ith entry being defined as By comparing (8) with (4), it is found that (8) can be regarded as a single sample received signal at a "new type vector sensor array" with the manifold matrixC Q .
Moreover, as discussed above, the manifold matrixC Q has two parts, in which the particle-velocity coarray manifold matrixC is determined by the response of the AVS and the spatial coarray manifold matrixQ is determined by the geometry of the array. The 16L 2 × 1 vectorr can be divided into 16 nonoverlapping subvectors, with each has a size L 2 × 1. In this way, each subvector would correspond to the same spatial coarray but to a different particle-velocity coarray component. For example, the ith subvector corresponds to the ith particle-velocity coarray component, i.e., the ith row ofC, and constitutes the ((i − 1)L 2 + 1)th to (i L 2 )th elements ofr. Let us denote the ith subvector ofr byr i , such that it can be extracted fromr and is expressed asr and Arrangingr i for i = 1, · · · , 16 in column-by-column format yields an L 2 × 16 augmented sample data matrix, which can be expressed in compact form as and Eq. (13) behaves as a multi-sample data augmented matrix "measured" by the spatial coarray. Also, the matrixR is constructed from the vectorization of the array covariance matrix by particle-velocity coarray augmentation. Here in this paper, the matrixR is called the PVCA matrix.

Identifiability of the PVCA Matrix
Using the subspace structure ofR, the signal subspace and the noise subspace matrices ofR can be obtained by performing the singular value decomposition onR and taking the corresponding left singular vectors. Therefore, the subspace-based algorithms, such as MUSIC [29] and ESPRIT [28], can be applied toR with the steering vector q(θ, φ) q * (θ, φ) ⊗q(θ, φ) for direction finding. In order to apply the subspacebased techniques, the signal part of the PVCA matrixR, i.e.,Q R sC T need to be of full-rank K . Indeed, the rank ofQ R sC T is dependent on the ranks ofQ, R s , andC individually. By properly designing the array topology, the matrixQ can be guaranteed to be full rank. Moreover, the diagonal matrix R s is always full rank. The problem remaining is to investigate the rank ofC, presented in the following theorem: The rank ofC is up to J = 13 when the AVS is deployed at or near a reflecting boundary and is up to J = 9 when the AVS is deployed in free space (far away from a reflecting boundary).
Proof See Appendix B.
By Theorem 1, the rank of the signal part ofQ R sC T is limited by min(rank(C),L), whereL denotes the maximum achievable rank ofQ. Since rank(C) ≤ min(J , K ), then In order for a subspace-based algorithm to work, it is required that The requirement (15) together with the constraint (14) indicates that the condition K ≤ J must be satisfied. Additionally, a general requirement for the subspace-based algorithms is that K ≤L − 1. These two constraints together lead to Therefore, by applying the subspace-based algorithms directly toR, up to 13 source signals can be identified for the reflecting boundary included scenarios, and up to 9 source signals can be identified for the reflecting boundary free scenarios.

Analysis of the PVCA
It is worth noting that in the PVCA, the matrixC is dependent only on the source directions, reflection coefficients, and the distance between the boundary and the AVS array, but not on the geometry of the spatial coarray. Hence, the PVCA is applicable to arbitrary array geometries and does not reduce the spatial coarray aperture. This flexibility constitutes a major advantage of the PVCA over the commonly used coarray processing techniques like the spatial coarray Toeplitz recovery technique [20], in which the coarray must be of a uniform structure and the spatial coarray aperture will be decreased after Toeplitz recovery. As the entire spatial coarray aperture can be exploited, the PVCA is applicable to deal with the underdetermined cases for identifying more sources than sensors. This can be readily justified because the spatial coarray consists of enough virtual sensors even for a "small" physical array of few sensors. For example, as depicted in Fig. 1, the spatial coarray of a 3 × 3 nonuniform planar array (in the x-y plane) would constitute a 7 × 7 "large" array. In this example, the PVCA can resolve up to 13 sources with the use of 9 physical sensors.
In establishing the results of Theorem 1, the reflection coefficient can take any values of R = ±1. This result is useful when considering some specific problems. For example, for the problem of direction finding with passive arrays mounted on the seabed [6], it is reasonable to model the interface between sea water and the packed sandy ocean bottom as the boundary between two fluids, one of which is absorptive [18]. In this way, according to [3], the reflection coefficient is given by where γ π/2−θ denotes the incident angle, η denotes the ratio between sand density and water density, and n denotes the refraction index. Moreover, the characteristic of absorption is modeled by assuming that the refraction index is complex-valued, i.e., n = n 0 (1 + jα), with α > 0. Note that, for a specific seabed interface, the reflection coefficients are dependent only on the source directions so that they would vary with the changing of the source directions, thereby ensuring the availability of the identifiability results in Theorem 1. The identifiability performance of the PVCA will be validated subsequently in Sect. 5 by considering the sandy ocean bottom scenario. For applications such as direction finding with a high frequency AVS array deployed near a vessel's hull [11], the reflecting surface can be modeled as the rigid boundary, in which the reflection coefficient takes a fixed value with R(θ ) = 1, ∀θ . In this case, there always exists extra linear dependence among the rows of the matrixC. Specifically, the normal velocity component is zero at d = 0, according to (1). Therefore, the number of the resolvable sources of the PVCA reduces to J = 6. Nonetheless, such reduction can be alleviated by setting the distance d = 0. It can be readily checked that for the case d = 0, the number of the resolvable sources of the PVCA is J = 10.
For applications such as direction finding with a low frequency AVS array deployed near a vessel's hull [11], the reflecting surface can be modeled as the pressure-release boundary, in which the reflection coefficient takes a fixed value with R(θ ) = −1, ∀θ . In this case, both the pressure and in-plane velocity components are zero at d = 0. This would cause rank(C) = 1, resulting in the failure of the PVCA. However, such failure can be overcome by setting the distance d = 0. Same to the rigid boundary case, the PVCA can identify J = 10 sources for the pressure-release boundary with d = 0.

Implementation Consideration
The derivation and analysis presented in this section are based on the covariance matrix R, whose exact value is generally unavailable in practice. To implement the PVCA for direction finding, the covariance matrix R may be replaced by its sample approximation as follows:R which is the maximum likelihood estimate of R. It is shown in [7] that, under the stationarity and ergodicity conditions,R will converge to R when the number of samples T tends to infinity. In summary, direction finding with PVCA is composed of the following main steps: (S1) Estimate the array covariance matrix R from the collected snapshots using (18). (S2) Perform vectorization on the estimate of R to get r. (S3) Perform row-exchanging using (8) to getr. (S4) ObtainR fromr using (13). (S5) Apply the MUSIC algorithm toR to obtain the azimuth-elevation angles of the sources.

Combination PVCA with Spatial Smoothing
Since the PVCA is built upon the exploitation of the coarray structure of the 4 × 1 AVS steering vector, there is a limitation on the number of the identifiable sources by applying the subspace-based algorithms directly to the PVCÃ matrixR. Nonetheless, this constraint can be readily mitigated by synergizing the PVCA with the spatial smoothing technique [4]. In order to apply the spatial smoothing technique, it is required that the spatial coarray is of a uniform structure, i.e., the spatial coarray manifold matrixQ is associated with a uniform array. This condition can be easily satisfied even if the physical array is of a nonuniform geometry. For example, it has been found that the spatial coarrays of the systemically designed nonuniform arrays, such as minimum redundancy arrays [22,26,40], coprime arrays [27,33,39], and nested arrays [8,13,24], exhibit the uniform structures. To focus only on the synergy between the PVCA and the spatial smoothing, the method for generating the spatial coarray with a uniform geometry from a systemically designed nonuniform array is not elaborated in this paper. Interested readers are referred to [8,13,22,24,26,27,33,39,40] for details. Now, letQ be the spatial coarray manifold matrix that is associated with aL ×L uniform square coarray in the x-y plane. To illustrate, an example of uniform square coarray withL = 7 is shown in Fig. 1. Obviously,Q is of sizeL 2 × K and is a subset (after some row exchanges) ofQ. To formulate the synergy between the PVCA and the spatial smoothing mathematically, the manifold matrixQ is represented as Q x Q y , whereQ x andQ y denote theL × K manifold matrices that are associated with the virtual sensors placed along the x-axis and y-axis, respectively. Based on the idea of the spatial smoothing technique, the entireL ×L spatial coarray is required to be partitioned into several overlapping subarrays, where each has the same size. These subarrays would thus have the translational invariant property from one another. Next, let the subarray size be L x × L y , where L x and L y represent the subarray dimensions along the x-axis and y-axis, respectively. Clearly, altogether (L −L x +1)×(L −L y +1) subarrays can be formed. With the preceding notations, the combination between the PVCA and the spatial smoothing is formulated in (19), whereL x L − L x + 1, L y L − L y + 1, J x denotes the x th to ( x + L x − 1)th rows of IL , and J y denotes the y th to ( y + L y −1)th rows of IL . Note that the assumption that all sensors lie in the x-y plane is not inherently necessary, but this will simplify the mathematical expressions in the analysis in order to obtain the identifiability results from (19).
It follows directly from (13) and (19) that the rank ofR satisfies rank(R) ≤ min(L xL y J ,L xL y K ).  Constraints (15) and (20) together lead to K ≤L xL y J . To apply the subspace-based techniques, it is required that K ≤L xL y − 1. Combining these two constraints yields where P =L xL y is the number of subarrays in performing spatial smoothing. Therefore, applying the subspace-based techniques directly toR, the maximum number of identifiable sources is K = min(L xL y − 1, P J ). With the combination to the spatial smoothing technique, the maximum number of identifiable sources offered by the PVCA can be multiplied by a factor P ≥ 2, which is a significant improvement over the results obtained for the use of the PVCA alone.

Computer Simulations
In the following simulations, a 4 × 4 nonuniform square AVS array, with sensors lying at the lattices generated by [0, 1, 4, 6]λ/2 in the x-y plane, is used as the measurement array. Both the reflecting boundary involved and the free space scenarios will be considered for evaluating the performance of the PVCA. The sandy ocean bottom boundary is considered for the reflecting boundary involved scenario. The corresponding reflection coefficient would follow the seabed model (17) with the typical values n 0 = 0.83, η = 0.27, and α = 0.1 [3]. The PVCA is compared to the vector field smoothing algorithm (labeled as "vector smoothing") [32] for the reflecting boundary involved scenario and to the tensor decomposition algorithm (labeled as "tensor decomposition") [34] for the free space scenario. In addition, the spatial coarray Toeplitz recovery technique (labeled as "Toeplitz recovery") [20] is modified to accommodate both the boundary involved and the free space scenarios and is included in comparisons. Unless otherwise stated, T = 1000 data samples are generated for implementing the algorithms.

Identifiability Performance
The identifiability performance of the PVCA is assessed in this example. First, consider a reflecting boundary involved scenario with the boundary lying in the plane z = λ/4. There are K = 13 sources, whose direction-cosines are listed in Table 2. The source signals are assumed to be equally powered with a 20 dB signal-to-noise ratio (SNR). By applying the MUSIC algorithm to the PVCA matrixR, the spatial spectrum shown in Fig. 2 is computed, where both 3-D mesh surface and the corresponding contour Table 3 Source Direction-Cosine Parameters (Free Space Scenario)  Table 3. The SNR remains at 20 dB.L x =L y = 3 is utilized to perform the synergy between the spatial smoothing and the PVCA to obtainR. The 3-D mesh surface and the corresponding contour plot of the spatial spectrum, which is computed by applying the MUSIC algorithm toR, are depicted in Fig. 3. One can observe that the synergy between the spatial smoothing and the PVCA successfully identifies these 18 sources. The identifiability results obtained from Figs. 2 and 3 are in consistent to our theoretical analysis presented in Sects. 3.2 and 4 .

Estimation Accuracy in Reflecting Boundary Involved Scenario
The estimation accuracy of the PVCA in the presence of a reflecting boundary is evaluated in this example. Two sources are located at θ 1 = 10 • , φ 1 = 20 • and θ 2 = 20 • , φ 2 = 30 • . The root-mean-squared error (RMSE) of source direction estimates is used as the metric for measuring the accuracy. By varying the SNR from 0 dB to 40 dB, the RMSE curves for θ 1 , θ 2 , φ 1 , φ 2 estimates of the vector smoothing, Toeplitz recovery, and PVCA algorithms are plotted in Fig. 4, where the corresponding CRB curves are also shown for the reference of performance benchmark. Here in Fig. 4, each RMSE datum is calculated from 500 independent trials and each CRB datum is computed by using Eqs. (42)-(51) in [12]. In general, for all the considered algorithms, the estimation RMSEs decrease linearly (in a semilog scale) with the increase in the SNR. Moreover, than the two comparative algorithms and the difference between the RMSE of the PVCA and the CRB is small.

Estimation Accuracy in Free Space Scenario
The estimation accuracy of the PVCA is compared to the tensor decomposition and the Toeplitz recovery algorithms for free space scenario in this example. The simulation parameter settings are the same as those in the last example. Figures 6 and 7

Conclusions
Using an AVS array, a new scheme has been proposed in this paper to estimate azimuth-elevation directions of multiple acoustic source signals. In this new scheme, the particle-velocity coarray augmentation (PVCA), which allows using the subspacebased techniques directly, is defined to increase the number of identifiable sources for There are two sources located at DOAs θ 1 = 10 • , φ 1 = 20 • and θ 2 = 20 • , φ 2 = 30 • . 500 independent trials are performed realizing underdetermined direction finding. In contrast to the widely used spatial coarray Toeplitz recovery technique, the PVCA processing does not reduce the efficient coarray aperture and is applicable to arbitrary array geometries. Moreover, the PVCA can be used in applications, where the array is deployed either near or far away from a reflecting boundary. For the reflecting boundary involved scenario, the PVCA can resolve up to 13 sources, while for the reflecting boundary excluded scenario, the PVCA can identify 9 sources at most. By applying to the recently developed nonuniform arrays like coprime arrays and nested arrays, the PVCA can be synergized with the spatial smoothing technique to get the number of resolvable sources multi-plied. Numerical simulations verify our theoretical results, thereby demonstrating the effectiveness of the proposed PVCA scheme.

Data Availability Statement
The data used to support the findings of this study are available from the corresponding author upon request.

Conflict of interest
We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, and there is no professional or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled, "Particle-Velocity Coarray Augmentation For Direction Finding with Acoustic Vector Sensors." Appendix A: Derivation of the Matrix P in (7) Define a P Q × P Q matrix U P×Q as where E pq is a P × Q matrix with all zeros except a 1 at the ( p, q)th position and F qp is a Q × P matrix with all zeros except a 1 at the (q, p)th position. Using the following three matrix prosperities given in [17] A we have and (I 4 ⊗ U 4×L 2 )(C * Q * ) ( Q C) = I 4 C * U 4×L 2 [( Q * Q) C] = (C * C) ( Q * Q).