DOA Estimation of Two-Dimensional Coherently Distributed Sources Based on Spatial Smoothing of Uniform Rectangular Arrays
Aiming at the direction-of-arrival (DOA) estimation of two-dimensional (2D) coherently distributed (CD) sources which are coherent with each other, we explore the propagator method based on spatial smoothing of a uniform rectangular array (URA). The rotational invariance relationships with respect to the nominal azimuth and nominal elevation are obtained under the small angular spreads assumption. A propagator operator is constructed through spatial smoothing of sample covariance matrices firstly. Then, combination of propagator and identical matrix is divided according to rotational operators, and the nominal angles can be obtained through eigendecomposition lastly. Realizing angle matching automatically, the proposed method can estimate multiple DOAs of 2D coherent CD sources without spectral peak searching and prior knowledge of deterministic angular signal distribution function. Simulations are conducted to verify the effectiveness of the proposed method.
In the field of array signal processing, traditional DOA estimation is based on point source models. In the real surroundings of radar and sonar systems, because of multipath propagation between receive arrays and targets, especially when the distances of targets and receive arrays are short, the spatial scattering of targets cannot be ignored, and the assumed condition of point source models is no longer valid. In such a condition, DOA estimation based on distributed source models has presented better accuracy . A distributed source can be regarded as an assembly of point sources within a spatial distribution where point sources can be called scatterers.
According to coherence of scatterers, distributed sources can be classified as coherently and incoherently distributed sources . Incoherently distributed (ID) sources are defined as scatterers within a source which are uncorrelated. On the contrary, coherently distributed (CD) sources are defined as scatterers within a source which are coherent. According to spatial distribution of sources, distributed sources can be classified as one-dimensional (1D) distributed sources and two-dimensional (2D) distributed sources. Under the assumption that scatterers of a target and receive array are not in the same plane, the 2D distributed sources should be more generally.
No matter CD or ID sources, most estimators suppose that signal from different sources is uncorrelated. In the field of underwater detection, the general process of detection is emitting narrowband pulse sound signal firstly, followed by receiving and analyzing the target’s backscatter signal. Thus, different CD sources are supposed to be coherent, which is more reasonable on account of the coherent multipath characteristics of underwater acoustic channel.
As to ID sources, utilizing different array configurations representative estimators have been proposed in [2–11]. Supposing that signals from different CD sources are coherent, sample covariance matrices are rank deficient. Thus, subspace-based algorithms [1, 2, 12–19] and ESPRIT class algorithms [20–23] which are based on eigendecomposition of sample covariance matrices cannot be applied for DOA estimation. PM class algorithms [24–27] based on linear operation of full rank sample covariance matrices are no longer applicable too. Considering signals from different CD sources are coherent, authors of  have proposed a DOA estimator by means of Toeplitz operation of sample covariance matrices. Authors of  have proposed a generalized MUSIC algorithm based on spatial smoothing. Both of the two algorithms above deal with 1D CD sources. Authors of  have proposed a general model of 2D CD sources coherent with each other which are defined as 2D coherent CD sources and proposed two spatial smoothing methods with double parallel linear arrays (DPLA).
In [1, 6–23], CD sources have been regarded as uncorrelated with each other; estimators cannot be applied for CD sources which are correlated. In , CD sources correlated with each other are discussed while sources are supposed to be 1D case. Literatures [10, 31] have considered 2D DOA estimators under uniform rectangular arrays, whereas they deal with ID sources and point sources, respectively. In , authors have considered 2D coherent CD sources and DOA estimators based on DPLA. In this paper, we are concerned on DOA estimation for 2D coherent CD sources under uniform rectangular array (URA). 2D coherent CD source and URA are introduced first. Then, the relationships of rotational invariance within and between subarrays have been derived under the assumption of small angular spreads. Afterwards, for decoherence of 2D coherent CD sources, spatial smoothing approach of URA is proposed, and a modified propagator is detailed for DOA estimation of 2D coherent CD sources. The method proposed need not angles matching and prior knowledge of deterministic angular signal distribution function. The simulation outcomes indicate that the proposed method is effective and advanced as for DOA estimation of 2D coherent CD sources.
2. Array Configuration and Signal Model
Figure 1 shows the URA configuration placed in the xoz plane. Arrays consist of M × K sensors separated by d meters along the direction of the x and z axis. Suppose that there are q 2D CD sources with nominal angles (θi, ϕi) () impinging the arrays, where θi is nominal azimuth angle and ϕi is nominal elevation angle of ith source, θi ∈ (0, π), ϕi ∈ (0, π). The noise is assumed to be additive Gaussian white with zero mean and uncorrelated between sensors.
Denote ymk(t) as the signal received by (m, k)th sensor, which can be expressed as follows:where is the impinging signal to the ith CD source, iL is the number of scatterers of the ith source, and (θil, ϕil) is the azimuth and elevation of the ilth scatterer of the ith source. Suppose different scatterers of the same source differ by one phase delay and a random amplitude gain. αil is complex gain of the ilth scatterer reflecting the reflection coefficient. E(αil) = αi. nmk(t) are noises received.
(θ, ϕ; ui) denotes the deterministic angular signal distribution function (ASDF) of the ith source which reflects the distribution of scatterers of a source and is related to the geometry and surface property of the source. A 2D ASDF is generally characterized by the parameter set ui = (θi, ϕi, σθi, σϕi) with four elements denoting the nominal azimuth, nominal elevation, azimuth angular spread, and elevation angular spread, respectively. Nominal azimuth and nominal elevation represent the center of the target which also can be expressed as DOA; azimuth spread and elevation spread indicating 2D spatial extension of a target are collectively called angular spreads.
For a Gaussian CD source, which means scatterers of the source obey Gaussian distribution, deterministic ASDF of the source can be expressed as
For a uniform CD source, which means scatterers of the source obey uniform distribution, deterministic ASDF of the source can be expressed as
Then, the received signal of (m, k)th sensor can be expressed as (see Appendix A)where si(t) = iLαi is the reflected signal of the ith CD source.
Apparently, when σθi = σϕi = 0, can expressed as a 2D Dirac function . Then, equation (4) can written aswhich is equivalent to the received signal model with regard to the point source. Most CD source models assume that scatterers within a source are coherent but scatterers between different sources are uncorrelated. According to characteristics of underwater acoustic channel, the assumption that scatterers between different sources are also coherent is more appropriate. So the reflected signal of the ith CD source can be written aswhere s1(t) reflects of time behavior of 1st source and ρi is the complex correlation coefficient between the ith source and 1st source. Define the generalized steering coefficient ηmk(θi, ϕi) of (m, k)th sensor as
Reflecting the responses of (m, k)th sensor to all CD sources, the generalized steering vector of (m, k)th sensor can be written as
Thus, the received signal of (m, k)th sensor can be written aswhere ρ = is the correlation coefficient vector. When the angular spreads of 2D CD sources are small and d/λ = 1/2, the following rotational invariance relationships can be obtained (see Appendix B):
3. Proposed Method
3.1. Subarray Signal Model
As shown in Figure 2, each subarray has M0 and K0 sensors along the direction of the x and z axes, respectively. Thus, subarrays smooth Ms = M − M0 + 1 times along the direction of the x axis and Ks = K − K0 + 1 times along the direction of the z axis.
Considering the (1, 1)th subarray, the receive vector of M0 sensors along the x axis can be expressed aswhere superscript 11 denotes the variable or vector with regard to (1, 1)th subarray. Define array as the closest and parallel to x axis of the (1, 1)th subarray. The array can be defined similarly as the (m, k)th subarray. is the noise vector of array .
Reflecting the responses of the array to all CD sources, is the generalized steering vector of , which can be defined as
According to equations (10)–(12), the received vector of closest kth and parallel to x axis array in the (1, 1)th subarray can be expressed aswhere Φz is the rotation operator, which can be written as
The received vector of the (1, 1)th subarray along the x axis can be expressed aswhere is the noise vector of the (1, 1)th subarray along x axis, which can be expressed as
Considering the (1, 1)th subarray and defining array as the closest and parallel to z axis of the (1, 1)th subarray, the receive vector of K0 sensors along the z axis can be expressed aswhere is the generalized steering vector of , which can be defined as
The received vector of the (1, 1)th subarray along the z axis can be expressed aswhere is the noise vector of the (1, 1)th subarray along z axis, which can be expressed as
Define as a vector selecting from 1 to M0K0 − M0 row of . is a vector selecting from M0 + 1 to M0K0 row of . and can be expressed aswhere (1:, M0K0 − M0:) denotes vectors from 1 to M0K0 − M0 row of . is noise vector selecting from 1 to M0K0 − M0 row of . is noise vector selecting from M0 + 1 to M0K0 row of . The generalized steering vector of can be expressed as
Considering the (m, k)th subarray, and are the sub vectors of , which can be expressed as follows:
Similarly, sub vectors and of can be obtained as
The generalized steering vector of can be expressed as
Considering the (m, k)th subarray, and are the sub vectors of , which can be expressed as
3.2. Spatial Smoothing and the Propagator Operator
Define the covariance matrices of the received signal of the (1, 1)th subarray aswhere , (·)H denotes the Hermitian transpose and N1, N2, N3, and N4 are covariance matrices of the noise vector. The spatial smoothing covariance matrices of the received signal can be obtained as follows:where
,, , and are spatial smoothing covariance matrices of the noise vector. The author of  has proved the spatial smoothing covariance matrices are full rank matrices when M0 > q, K0 > q and MS ≥ q, KS ≥ q.
Construct matrix B aswhere B1 is a q × q dimensional nonsingular matrix and B2 is a (4M0K0 − 4M0 − q) × q matrix. There exists a q × (4M0K0 − 4M0 − q) dimensional propagator operator P satisfying .
Letwhere Iq×q is q × q identity matrix. B can be expressed as follows:
Sample covariance matrices with N snapshots of the (m, k)th subarray can be defined as
Thus, spatial smoothing of sample covariance matrices can be obtained as
Combine ,,, and into
Divide (4M0K0 − 4M0) × (M0K0 − K0) dimensional matrix into q × (M0K0 − K0) dimensional matrix G and (4M0K0 − 4M0 − q) × (M0K0 − K0) matrix H:
The propagator operator P can be obtained aswhere (·)+ denotes the pseudoinverse. Divide E into four (M0K0 − M0) × q dimensional matrices E1, E2, E3, and E4.
From this equation, the following relationship can be obtained:
So, we have
We can obtain eigenvalue μi () and it’s corresponding eigenvectors ξi of Ω1 by means of eigendecomposition. From equations (49) and (50), we can conclude that the corresponding eigenvalue of Ω1 and Ω2 have the same eigenvector; so the eigenvector of Ω2 can be obtained as follows:where (./) denotes element-wise division operation. When d/λ = 1/2, nominal elevation angle ϕi and nominal azimuth angle θi of ith source can be obtained as follows:where denotes argument of complex variable.
3.3. Computational Procedure Complexity Analysis
Now, our algorithm can be summarized as follows: Step 1: compute spatial smoothing of sample covariance matrices using equations (39)–(43). Step 2: estimate the propagator operator P using equation (46). Construct E and divide E into E1, E2, E3, E4 from equation (47). Step 3: find eigenvalue μi () and it’s corresponding eigenvectors ξi through eigendecomposition of Ω1. Step 4: obtain eigenvalue νi of from equation (50). Step 5: calculate the nominal azimuth θi and nominal elevation ϕi from equation (51).
The computation cost of the proposed method mainly consists of three parts: the calculation of the sample covariance matrix which is the number of snapshots times the number of spatial smoothing multiply multiplication number with regard to calculation of covariance matrices of subarray, calculation of Ω1 and Ω2 where computation cost mainly lies in pseudoinverse operation O[(2M0K0 − 2M0)3], and eigendecomposition of Ω1 and Ω2O(q3). Suppose 2D CD sources are incoherent, which means decoherence by spatial smoothing is unnecessary. Then, computation cost mainly consists of calculation of the sample covariance matrix , pseudoinverse operation O[(2MK − 2M)3], and eigendecomposition O(q3). It can be concluded that using spatial smoothing to realize decoherence of 2D coherent CD sources does not change the magnitude of computational complexity in essence.
4. Simulation Results
In this section, five simulation experiments are conducted to verify the effectiveness of the algorithm we proposed. All simulation experiments are based on array configuration as shown in Figure 1. The distance of adjacent sensor is set at .
denotes root mean squared error (RMSE) of the ith source, which can be expressed aswhere Mc is the Monte Carlo simulations number which is set at 100. and are the estimated nominal azimuth and nominal elevation of ith source in ςth Monte Carlo simulation.
In the first example, we investigate the performance of the proposed method versus three 2D CD sources which are uncorrelated with each other. The first source is Gaussian CD source with parameter set (30°, 45°, 2°, 2°). The second and third sources are uniform with parameter sets (50°, 45°, 2°, 2°) and (50°, 60°, 2°, 2°). The number of snapshots is set at 200. The subarray number parameter is Ms = Ks = 4, subarray elements parameter is M0 = K0 = 4. RMSE is the mean of three sources. As shown in Figure 3, the proposed method presents better estimation as the SNR increases. Figure 3 also shows estimation of the DSPE  which uses L-shaped arrays containing sensors in x and z axes, and the ESPRIT  method uses DPLA. It can be concluded that DOA estimation of all the three methods is effective with respect to 2D CD sources which are uncorrelated with each other.
In the second example, we investigate the estimation of 2D coherent CD sources versus SNR and number of snapshots. We consider three full coherence sources with equal power. The first source is Gaussian with parameter set [30°, 45°, 2°, 2°]. The second and third sources are uniform with parameter sets [50°, 45°, 2°, 2°] and [50°, 60°, 2°, 2°]. Ms = Ks = 4 and M0 = K0 = 4. RMSE takes mean values of the three sources. Figure 4(a) shows RMSE with SNR varying from 0 dB to 30 dB while the number of snapshots is set at 200. Figure 4(b) shows RMSE with the number of snapshots ranging from 100 to 1000 while SNR is fixed at 15 dB. It can be observed that the proposed algorithm performs better as SNR or number of snapshots increases. Nevertheless, the DSPE  and ESPRIT  present big errors as shown in Figures 4(a) and 4(b). It can be concluded that utilizing traditional methods for 2D coherent CD sources is invalid.
In the third example, we investigate the performance of the proposed method versus angular spreads. We consider two scenarios. One scenario has three full coherence Gaussian sources with equal power, and the nominal angles are (30°, 45°), (50°, 45°), and (50°, 60°). The other has three uniform sources with the same nominal angles as the first scenario. To simplify the analysis, we assume that azimuth angular spread is equal to elevation angular spread within each source, and angular spreads of three sources are the same; σ is used to replace and for convenience. The number of snapshots is set at 200 and SNR is 15 dB. Ms = Ks = 4 and M0 = K0 = 4. RMSE of the two trails is defined as the mean values of three sources. From Figure 5, it can be observed that estimation errors increase as angular spreads increase with respect to both Gaussian and uniform sources. RMSE reaches 0.027° as angular spreads increase to 5° and reaches 0.4° as angular spreads increase to 10°, which still is a satisfactory result. It can be seen that the proposed method has a good performance under the small angular spreads condition. It is noteworthy that when the angular spreads is 0°, 2D CD sources degenerate into point sources; this experiment also manifests that the method proposed can also be applied to point sources.
In the fourth example, we investigate the performance of the proposed method versus subarray parameters. We consider three full coherence sources with equal power. The source parameter sets are the same as the second example. The number of snapshots is set at 200, and SNR is 15 dB. For the convenience of analysis, Ms is supposed equal to Ks, and M0 is supposed equal to K0. Figure 6 displays the performance of the proposed method versus the subarray number parameter Ms with M0 = 4, while Figure 7 shows the performance versus subarray elements parameter M0 with Ms = 4. From Figures 6 and 7, it can be concluded that the proposed method can estimate DOA of the 2D coherent CD sources effectively when M0 > q, K0 > q, MS ≥ q, and KS ≥ q and provides better performance as subarray number or subarray elements increases, but RMSE changes slightly when M0 or Ms increases to a certain extent.
In the fifth example, we investigate the performance of the proposed method near the boundary region in comparison with the method  utilizing DPLA. We consider double full coherence Gaussian sources with equal power. The number of snapshots is set at 200, and SNR is 15 dB. Angular spreads are set at 2°. In the first trail, nominal elevations of two sources are set at 30° and 50°, and nominal azimuths of two sources are varying from 0° to 180° simultaneously. Similarly, in the second trail, nominal azimuths of two sources are set at 30° and 50°, and nominal elevations of two sources are varying from 0° to 180° simultaneously. Ms = Ks = 3 and M0 = K0 = 3, and the total sensor number of URA is 25. DPLA is also set with total sensor number as 25. RMSE takes the mean values of the two sources. Figures 8(a) and 8(b) show the first and second trails, respectively. As shown in Figures 8(a) and 8(b), the method proposed has a better performance near the boundary regions and more significant near the elevation boundary region. Aperture sizes of URA are balanced in azimuth and elevation dimensions. Compared with URA, DPLA has smaller aperture in the elevation dimension. Therefore, errors of elevation estimation near the elevation boundary region are greater; consequently, estimation deterioration will be more significant in these areas as solving elevation and azimuth of both methods are coupled.
In this study, we have considered the problem of estimation of 2D coherent CD sources utilizing URA. The rotation invariance relationships within and between subarrays have been deduced under the small angular spreads assumption. Decoherence can be realized by virtue of 2D spatial smoothing of URA. Propagator method of sample covariance matrices based on spatial smoothing has been introduced in detail. Simulation investigates experiment conditions containing SNR, number of snapshots, angular spreads, and subarray parameters. Simulation outcomes indicated that the proposed method is effective for DOA estimation of 2D coherent CD sources and has a better performance near the boundary regions. The method proposed is also effective for the 2D CD sources which are incoherent as well as point sources.
A. Integral Expression of Received Signal
Considering equation (1), angle (θil, ϕil) which is DOA of ilth scatterer of the ith source can be regarded as 2D random variables, and when iL is a large number, converges to mean value of according to the law of large numbers. Then, we have
As the scatterers of the ith source are subject to the density distribution ,
Then, ymk(t) can be expressed as
B. Rotational Invariance Relationships
Change the variables () for θ and () for ϕ, where and are the small deviations of θi and ϕi. Thus, cos θ sin ϕ and cos ϕ can be approximated by the first term in the Taylor series expansions. Consider the relation between and :
Because of the following relationshipwe have
Consider the relation between and :
Noticing the following relationship
So, we have
Similarly, we have
The data used to support the findings of this study are available from email addresses [email protected].
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Tao Wu and Yiwen Li contributed equally to this work. Tao Wu and Yiwen Li designed the algorithm scheme. Tao Wu and Xiaofeng Zhang designed the software and performed the experiments. Tao Wu and Chaoqi Fu wrote the first draft. Yijie Huang and Qingyue Gu performed proofreading and editing. All authors read and approved the final manuscript.
This research was funded by the National Natural Science Foundation of China (Grant numbers 61471299 and 51776222).
P. F. Ma and S. G. Yan, “Estimation of the DOA of coherent distributed sources,” Technical Acoustics, vol. 33, pp. 817–821, 2007.View at: Google Scholar