Abstract

Nested arrays are sparse arrays composed of subarrays with nonuniform sensor spacing. Compared with traditional uniform arrays, nested arrays have more degree of freedoms (DOFs) and larger apertures. In this paper, a nested array has been proposed as well as a direction-of-arrival (DOA) estimation method for two-dimensional (2D) incoherently distributed (ID) sources. A virtual array is firstly obtained through vectorization of the cross-correlation matrix of subarrays. Sensor positions of the virtual array and the optimal configuration of the nested array are derived next. Then rotational invariance relationship for generalized steering matrix of the virtual array with respect to nominal azimuth is deduced. According to the rotational invariance relationship, sparse representation model under l1-norm constraint is established, which is resolved by transferring the objective function to second-order cone constraints and combining a estimation residual error constraint for receive vector of the virtual array. Simulations are conducted to investigate the effectiveness of the proposed method in underdetermined situation and examine different experiment factors including SNR, snapshots, and angular spreads as well as sensor number of subarrays. Results show that the proposed method has better performance than uniform parallel arrays with the same number of sensors.

1. Introduction

In the field of array signal processing, traditional DOA estimation of targets is based on point source models. Point source models suppose signal propagation from a target to receive array is one straight path; geometry of the target can be ignored, which can be considered as a point. In the surrounding of underwater acoustics detection, on one hand, there exist reflection paths from seabed and sea surface between a target and receive array; on the other hand, when the distance from the target to receive array is short, as many parts of the target reflect signal and spatial distribution and the geometry of the target cannot be ignored; signal propagations from the target to receive array cannot be supposed to a single path. Point source models cannot describe spatial properties of sources under such circumstance; estimations based on point source models are no longer effective, which should be replaced by distributed sources [1]. Distributed sources can be regarded as an assembly of point sources which can be called scatterers within a spatial distribution.

A distributed source may be coherently distributed (CD) source or incoherently distributed (ID) source according to the coherence of scatterers of a target. A coherently distributed (CD) source means that scatterers of the target are coherent. Spatial distribution of scatterers of a CD source is described by deterministic angular signal distribution function (ASDF). On the contrary, an incoherently distributed (ID) source means that scatterers of the target are uncorrelated. Spatial distribution of scatterers of an ID source is described by angular power density function (APDF). Both ASDF and APDF can be generally modeled as Gaussian or uniform distribution according to the geometry and spatial distribution of a target. In this paper, ID sources are considered.

Distributed sources may be one-dimensional (1D) or two-dimensional (2D) which depends on whether impinging signals of scatterers and arrays are in the same plane. As to a 1D source, parameters of ASDF or APDF include nominal angle and angular spread. Nominal angle can be called nominal DOA representing the center of a target while angular spread represents spatial extension of the target. As to a 2D source, parameters contain nominal azimuth, nominal elevation, azimuth spread, and elevation spread. Nominal azimuth and nominal elevation represent center of a target which also can be expressed as DOA; azimuth spread and elevation spread indicating 2D spatial extension of the target are collectively called angular spreads. As impinging signals of scatterers and arrays are not in the same plane practically, 2D sources are more general and reasonable.

As to CD sources, utilizing different array configurations representative estimators have been proposed in [210], which are mostly based on rotational invariance relationship derived by Taylor approximation of generalized steering vectors. For ID sources, several algorithms have been presented and can be roughly divided into four categories: subspace-based algorithms developed from multiple signal classification (MUSIC) such as distributed signal parameter estimator (DSPE) [1] and dispersed signal parametric estimation (DSPARE) [11], generalizations of Capon’s methods [1214] involving high-order matrix inversions and spectral searches, the maximum likelihood (ML) approaches [1517] requiring high computational complexity but with better accuracy, and covariance matching estimation techniques (COMET) [1821] with lower computational complexity than ML approach but the same large sample behavior.

The aforementioned estimations of ID sources suppose the sources and receive arrays are in the same plane, which are special cases modeled as 1D ID sources. Involving more parameters, there have been relatively few studies on estimation of 2D ID sources. An extension of COMET has been proposed in [22] which employs alternating projection principle separating source powers and noise variances from objective functions and then estimates parameters by a four-dimensional nonlinear optimization. Utilizing double parallel uniform linear arrays, a TLS-ESPRIT-like approach has been proposed for DOA estimation in [23], where nominal elevations are firstly estimated via rotational invariance relationship derived by the first-order Taylor series expansions of the generalized steering vectors; then nominal azimuths are obtained through 1D searching. An estimator based on uniform rectangular arrays has been proposed in [24], where DOAs are estimated through an ESPRIT like algorithm by virtue of rotational invariance relationship of receive vectors which are described by generalized signal vector and generalized steering vector composed of nominal steering vector and its first-order partial derivatives.

Most distributed sources estimation techniques mentioned above are based on uniform arrays such as uniform linear arrays, uniform rectangular arrays, or circular arrays, which require the distances between sensors less than or equal to half of impinging signal wave. For achieving better accuracy and resolved more sources, estimation should be performed with larger aperture, which needs to increase number of sensors. Consequently, this would lead to an increase in complexity and cost of sonar system. The recent proposed nested arrays [25] possess larger apertures and more degree of freedoms (DOFs) compared with uniform arrays with the same number of sensors. DOA estimations for point sources under different kinds of nested arrays have been proposed in [2630]. As for distributed sources, based on a linear nested array, a spatial smoothing with spectral search technique has been proposed in [8] for 1D distributed sources supposing a priori knowledge of the angular spreads is known, which is impractical in practice. Using the same array configuration, the authors of [21] have proposed annihilating filter technique combined with structured low rank approximation for 1D distributed sources without a priori knowledge of the angular spreads, where the aperture of virtual array has not been fully utilized as half of the sensors for estimation of DOA and the other half for angular spreads.

In this paper, DOA estimation for 2D ID sources is proposed based on a 2D nested array composed of double parallel uniform linear subarrays. According to the concept of difference coarray, the cross-correlation matrix of subarrays is firstly vectorized by Khatri-Rao product to obtain a hole-free virtual array; and the closed form expression for sensor positions of the virtual array is presented. Then rotational invariance relationship with regard to nominal azimuth is derived under the assumption of small angular spreads. Sparse representation with l1-norm constraint for DOA estimation of 2D ID sources is established next. Objection function is then converted to second-order cone constraints to avoid determining the regularization parameters. Based on estimation residual error for receive vector of virtual array, a robust constraint is derived. Followed by the fact that nominal azimuths are resolved by the sparse representation framework, paired nominal elevations are obtained by least-square techniques. The proposed method can detect 2D ID sources more than sensors without angles matching and without information of APDF of sources. Simulations are conducted to investigate the effectiveness and superior of the proposed method lastly.

2. The Nested Array and Source Model

Generally, nested arrays are composed of several levels of uniform linear subarrays; the distances between sensors within the first level subarray are less than or equal to half of impinging signal wave; distances between sensors within the rest subarrays are far greater than the first level subarray. As shown in Figure 1, the proposed nested array is composed of double parallel uniform subarrays on xoz plane X1 and X2. Subarray X1 is parallel to x axis and consists of the middle sensor placed on z axis and M1 sensors located on both sides. The sensor number of subarray X1 is N1=2M1+1; the distance between sensors of subarray X1 is d. Subarray X2 is on x axis with M2 sensors located on both positive and negative semiaxes. The sensor number of subarray X2 is N2=2M2; the distance between sensors of subarray X2 is D=N1d. The distance between subarrays X1 and X2 is also d. Suppose that there are q 2D ID narrowband sources with nominal angles (θi, φi) (i=1,2,⋯,q) impinging into the nested array. θi denotes a spatial nominal azimuth between x axis and propagation path of the ith source; φi is nominal elevation of the ith source. θi, φi. λ is the wavelength of the impinging signal.

Sensor positions set of subarrays X1 and X2 can be expressed, respectively, as

N1×1 dimensional receive vector of subarrays X1 and N2×1 dimensional receive vector of subarray X2 can be expressed aswhere and are additive noise vectors which are assumed to be uncorrelated between sensors and Gaussian white with zero mean. a1(θ,φ) and a2(θ,φ) denote the steering vector of subarrays X1 and X2 with respect to point source. si(θ,φ,t) is the complex random angular signal density of the ith distributed source representing the scatterer intensity of the source from angle (θ, φ) at the snapshot index t. Unlike a point source, signal of a distributed source exists not only in the single direction (θi, φi) but also in a spatial distribution around (θi, φi). An ID source means different scatterers from one target generate uncorrelated multipath signals. Therefore, one direction from si(θ,φ,t) is uncorrelated with other directions from si(θ,φ,t); we havewhere is power of the ith source, δ(•) is the Kronecker delta function, (•) means complex conjugate, and f(θ,φ; u) is the source’s normalized APDF which reflects the distribution of scatterers of a source and is related to the geometry and surface property of the source. APDF can be modeled as 2D Gaussian, uniform, or any other distribution function according to different characteristic of an ID source, which is determined by parameter set ui=[,,,] denoting nominal azimuth, nominal elevation, azimuth spread, elevation spread, respectively, satisfying

For a Gaussian ID source, APDF can be expressed as

For a uniform ID source, APDF can be expressed as

3. Proposed Method

This section consists of five parts. First of all, a virtual array is obtained by difference of sensor positions in two subarrays. Then, optimal configuration of the nested array is analyzed and the rotational invariance relationship of generalized steering matrix of the virtual array is derived. In the next part, DOA estimation via sparse representation is established and solution process is detailed. Lastly computational procedure is summarized and complexity analysis is analyzed with comparison of two existing methods for 2D ID sources.

3.1. Description of Virtual Array

The main advantage of nested array lies in increasing DOFs by forming virtual array through subarrays. According to the assumption of ID sources, different sources are uncorrelated. The N1×N2 dimensional cross-correlation matrix of subarrays X1 and X2 can be expressed as follows:

Vectorizing R by Khatri-Rao product, we obtainwhere rec(•) denotes vectorization of a matrix, denotes Khatri-Rao product, and (•)H denotes Hermitian transpose. A virtual array can be obtained by vectorization of R which is shown in Figure 2. Sensor positions of the virtual array are determined by difference of sensor positions in subarrays X1 and X2. The set of sensor positions of the virtual array can be expressed as

From (1) and (2) we obtain sensor position set of the virtual array aswhich indicates that the virtual array has 2(2M1M2+ M2) virtual sensors without repetitive positions.

It is noteworthy that elements in are not arranged in accordance with sensor order in virtual array. Suppose r is receive signal vector from left side of virtual sensor located on [(0.5-2M1M2-M2)d, 0, d] to the right side of sensor located on [(2M1M2+M2-0.5)d, 0, d]. Define a block diagonal matrix J for sorting order of , which is composed of N2 number of blocks which are N1 dimensional inverse identity matrix. J can be expressed as follows:

Then, r can be written aswhere b(θ,φ) denotes the steering vector of subarray virtual array with respect to point source, which can be written as

3.2. Optimal Configuration of the Nested Array

From (2), DOF of the virtual array is 2(2M1M2+M2); sensor number of the nested array is N=2(M1+M2)+1. It can be concluded that DOF of the proposed nested array can achieve O(N2) which is far more than number of sensors. We investigated the most DOF can be achieved under the given sensor number N; the optimization problem can be expressed as follows:

Equation (16) can be solved through minimization of Lagrange function as follows:

Taking the derivative of (17) with regard to M1, M2 setting the result equal to 0, we obtain

As M1, M2, and (N-1)/2 must be integers, optimal solutions are all fractions and need to be modified, which are listed in Table 1.

3.3. Rotational Invariance Relationship of Virtual Array

In the point source models case, steering vectors are used to describe response of an array. Nevertheless, as for distribute sources, generalized steering vectors or generalized steering matrices [5, 23] represent response of arrays. Define as the generalized steering vector of the virtual array with respect to the ith distributed source, which can be expressed aswhere the kth element of means the response of kth sensor to the ith distributed source. Generalized steering matrix of the virtual array can be expressed as

r can be expressed as where is the power vector with following expression:

Suppose the angular spreads are small and d/λ is set at 1/2; we have the following relationship (see Appendix A)where [•]k represent the kth element of vector. So we can obtain the following rotational invariance relationship:where is the rotation matrix, Bk is the kth row of B representing the response of kth sensor with respect to distributed sources.

3.4. DOA Estimation via Sparse Representation

Representing the response of first sensor to distributed sources, B1 is the first row of B. Change variables for θ and for φ, where and are the small deviations of and . B1 can be regarded as Hadamard product of three vectors which are nominal azimuth coefficient vector, nominal elevation coefficient vector (φ), and real coefficient vector g. B1 can be expressed as (see Appendix B)where ⊕ represent Hadamard product. (φ) and g can be written as where . According to the rotational invariance relationship (24), the kth row of can be obtained as

Define a nominal azimuth coefficient matrix as where

Thus, receive vector of virtual array can be expressed aswhere

Divide the nominal azimuth space into Q partitions. Suppose set = contains all the possible nominal azimuths. Sparse representation of receive vector r can be expressed aswhere is sparse coefficient vector. The jth element of vector equals the ith element of P if , (j=1,,Q; i=1,,q). Therefore, can be obtained from the positions of nonzero elements of . Though l0 norm is a direct measurement of sparsity of a vector, l0 norm constraint optimization is a NP-hard problem; we use l1 norm constraint to replace l0 norm constraint. Thus, DOA estimation via sparse representation under l1 norm constraint can be expressed as

The cross-correlation matrix R can be replaced by sample cross-correlation of subarrays with L snapshots as follows:

Then we have the vectorization of receive signal vector of virtual array as follows:where is estimation residual error of virtual array, which has the following probability distribution (see Appendix C) where AsN(μ, C) represents multidimensional normal distribution with mean value μ and covariance matrix C. ⊗ is Kronecker product. R1 and R2 are covariance matrices the subarrays X1 and X2, respectively. According to the nature of normal distribution, we have where W=, Asχ2(4M1M2+2M2) represents Chi square distribution with DOF as (4M1M2+2M2). Then we introduce a parameter η to restrain Δr as follows:

Let (40) hold with a high probability p=0.999. η can be resolved from chi2inv(p, 4M1M2+2M2) function of Matlab. Then (34) can be converted into the following robust formwhere 1 is Q×1 dimensional vector whose elements all are 1, (41) can be resolved by CVX. R1 and R2 can be replaced by sample covariance matrix with L snapshots.

Through tracking nonzero position of coefficient , θi can be resolved. Next we can obtain the estimated coefficient vector P as , where (•)+ denotes pseudo inverse operation. Thus, nominal elevation of the ith source can be obtained as follows:

3.5. Computational Procedure and Complexity Analysis

Now, our algorithm can be summarized as follows.

Step 1. Compute sample cross-correlation matrix R and sample covariance matrices R1 and R2 using (35) and (42).

Step 2. Vectorize R from (10) and reorder elements of the vector to obtain according to (36).

Step 3. Calculate covariance matrix of residual error vector according to (37) and determine upper bound of constraint of η.

Step 4. Discretize nominal azimuth and calculate nominal elevation θi from (41).

Step 5. Calculate nominal elevation φi from (43).

We analyze the computational complexity of the proposed method in comparison with COMET [22] which can be applied for uniform 2D arrays and Zhou’s algorithm [23] using double parallel arrays for 2D ID sources. For comparison, the two methods are supposed to be applied to double parallel arrays with the same sensor number as the proposed nested array. COMET [22] is a method of four-dimensional nonlinear optimization under alternating projection algorithm framework, its computational cost is mostly composed of the calculation of the sample covariance matrix which needs O(LN2) and the alternating projection technique with respect to cost functions which is O(3N3q3+ 3N2). Zhou’s algorithm calculates nominal elevation by a TLS-ESPRIT algorithm and finds nominal azimuth by 1D searching. Computational cost of Zhou’s algorithm [23] mainly contains calculation of the sample covariance matrix O(LN2), eigendecomposition, and inversion of the sample covariance matrix O(2N3) and 1D searching O(N3). Computational cost of the proposed algorithm mainly contains calculation of the sample cross-correlation matrix and sample covariance matrix O(LN2), solving second-order core programming O[(4M1M2+2M2)3Q3] and solving the nominal elevation by pseudo inversion O(N3). With more DOF the computational cost of the proposed algorithm is higher than two methods.

4. Results and Discussion

In this section, the effectiveness of the proposed algorithm is investigated through four numerical simulation experiments. We set the nested array with M1=2 and M2=3, the total sensor number is N=11, d=λ/2. The nominal azimuth space is divided from 0° to 180° in step of 1°. Root mean squared error (RMSE) is applied for the evaluation of the performance of estimation. The RMSE of the ith source with respect to DOA is defined aswhere and are the estimated nominal azimuth and estimated nominal elevation of the ith ID source, respectively. The superscript denotes the estimated parameters in th Monte Carlo run. Mc is the number of Monte Carlo simulations. θi and φi are the true nominal azimuth and nominal elevation, respectively.

In the first experiment, we examine the performance considering underdetermined scenario. 2D ID sources to be estimated are distributed from (15°, 15°) to (165°, 165°) uniformly. The center of kth sources is located on [15°+(k-1)10°, 15°+(k-1)10°]. APDF of sources are Gaussian when k is even while sources are uniform when k is odd. Angular spreads are all set at 2°, SNR is set at 5db, and number of snapshots is 200. 100 independent Monte Carlo simulations are run to get the result. Figure 3 shows the positions of given 16 sources as well as average estimated positions of the 16 sources, which indicate that the proposed method can estimate sources with different APDF effectively in underdetermined scenario.

In the second experiment, we examine the performance of the proposed method versus angular spreads of ID sources. As the rotational invariance relationship of generalize steering matrix is obtained under the small angular spreads assumption, it is necessary to investigate the influence of angular spreads on the estimation of the proposed algorithm. We consider both Gaussian and uniform sources, respectively. SNR is set at 5dB, number of snapshots is 200, and Mc=100. Sources are centered at (30°, 45°) and (150°, 70°) with both angular spreads ranging from 0° to 10°. RMSE takes averages of two sources. Figure 4 shows that RMSE reaches 0.6 as angular spreads are 5° and RMSE reaches 0.8 when angular spreads are 10° for Gaussian sources, while RMSE reaches 0.75 and 1.1 when angular spreads reach 5° and 10° for uniform sources, respectively. It can be concluded that estimation accuracy deteriorates for both Gaussian and uniform sources as the angular spreads increase. The estimated results are still satisfactory within 5° and the proposed method shows robustness under small angular spreads.

In the third experiment, we examine the performance of the proposed method versus SNR and number of snapshots. Two Gaussian ID sources with parameters [30°, 45°, 2°, 2°] and [150°, 70°, 2°, 2°] are set to be estimated. =100. Figure 5(a) shows RMSE of estimation with SNR ranging from -5dB to 30dB while number of snapshots is set at 200. RMSE takes averages of and of two sources. Figure 5(b) shows the influence of number of snapshots with SNR set at 5 dB. Sources are also estimated by COMET [22] and Zhou’s algorithm [23] based on two kinds of double parallel uniform linear arrays for comparison. One is composed of 12 sensors with 6 sensors separated by λ/2 meters in each subarray, which has one more sensor than nested array we proposed. The other has 30 sensors with each subarray containing 15 sensors separated by λ/2 meters. Physical aperture of the parallel uniform linear array with 30 sensors is similar to the proposed nested array. As can be seen from Figure 5, estimation accuracy increases as SNR or number of snapshots increases. As the aperture of nested array extended by the virtual array, the proposed method shows better performance than the two methods with the same number of sensors. Considering the parallel uniform linear array with 30 sensors, as shown in Figure 5(b), RMSE of COMET method is smaller than the proposed method when number of snapshots exceeds 200 and bigger than the proposed method when number of snapshots is less than 100. COMET separates unknown variables of each source based on alternating projection technique and formulates equations set of unknown variables from covariance matrix of receive vectors. Using more sensors, COMET shows better performance if the samples are sufficient. In the case of number of snapshots less than 100, the proposed method performs better as applying a robust sparse representation. In summary, with the same aperture, the proposed method has a better performance than uniform arrays in less snapshots circumstance; with same number of sensors, the proposed method has a better performance than uniform arrays under the same experimental conditions.

In the fourth experiment, we investigate the impact of sensors number of subarrays on estimation. Two Gaussian ID sources with parameters [30°, 45°, 2°, 2°] and [150°, 70°, 2°, 2°] are set to be estimated. Mc=100, SNR=5dB, number of snapshots is 200. Figure 6 shows four curves representing 4 sensor number settings. The distances between adjacent sensors d are all set at λ/2. Estimation accuracy will be improved as sensor number in subarrays increases; but this improvement is not significant when number of sensors in subarrays increases to a certain extent. Rotational invariance relationship is an approximate resolution under small angular spreads assumption; when sensors number is bigger, though aperture is larger, expression of receive signal with respect to later sensors in virtual array is a result of more transitive relationships, which means existing more errors.

5. Conclusions

In this paper a 2D nested array is proposed aiming at DOA estimation for 2D ID sources. Through vectorization of cross-correlation matrix by Khatri-Rao product, a virtual array is obtained. The close form of sensor position of virtual array and the optimal configuration of the nested array are derived. Rotational invariance relationship of generalized steering matrix with regard to nominal azimuth is deduced under the small angular spreads assumption. Sparse representation of DOA estimation for 2D ID sources based on the nested array is established and the solution process is detailed. Simulations investigate the influence of experiment conditions, angular spread, and sensor number on estimation. Results indicate that the proposed method can estimate more sources than sensors without information of APDF of sources and has better performance than uniform parallel arrays with the same sensor number.

Appendix

A. Rotation Invariant Relationship

Check []k and []k+1 of (23). Replace variables with θ and with φ, where and are the small deviations of and . Thus, and can be approximated by the first term in the Taylor series expansions. []k and can be, respectively, expressed as follows:

Under the assumption of small angular spreads, it follows that , so []k+1 can be written as

Then (23) can be obtained.

B. Hadamard Product Expression of B1

Considering [B1]i of (25), which is the ith element of the first row of B, according to (A.1), [B1]i can be expressed as follows:where [(φ)]i is the ith element of (φ) described by (26). Thus, (25) can be obtained.

C. Covariance Matrix of Estimation Residual Error

Denote (N1N2×1) dimensional as the estimated vector of which is vectorization of sample cross-correlation matrix of subarrays. The ith (N1×1) subvector of is given by

The signals are supposed independent from snapshot to snapshot; the following relation is easily obtainedwhere [•]ji denotes the element of jth row and ith column within the matrix. Equation (C.2) results in

Then, we have

Data Availability

The data used to support the findings of this study are available from email addresses [email protected] and [email protected].

Conflicts of Interest

The authors declare no conflict of interest.

Authors’ Contributions

Wu Tao and Yiwen Li designed the algorithm scheme. Wu Tao and Zhengxin Li wrote software and performed the experiments. Wu Tao wrote the first draft. Yijie Huang and Jiwei Xu made proofreading and editing. All authors read and approved the final manuscript. Tao Wu and Yiwen Li contributed equally to this work.

Acknowledgments

This research was funded by the National Natural Science Foundation of China grant number 61471299 and 51776222.