#### Abstract

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.

#### 1. Introduction

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 [1]. 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 [1]. 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 [28] have proposed a DOA estimator by means of Toeplitz operation of sample covariance matrices. Authors of [29] have proposed a generalized MUSIC algorithm based on spatial smoothing. Both of the two algorithms above deal with 1D CD sources. Authors of [30] 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 [24], 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 [30], 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 *i*th source, *θ*_{i} ∈ (0, π), *ϕ*_{i} ∈ (0, π). The noise is assumed to be additive Gaussian white with zero mean and uncorrelated between sensors.

Denote *y*_{mk}(*t*) as the signal received by (*m*, *k*)th sensor, which can be expressed as follows:where is the impinging signal to the *i*th CD source, *iL* is the number of scatterers of the *i*th source, and (*θ*_{il}, *ϕ*_{il}) is the azimuth and elevation of the *il*th scatterer of the *i*th source. Suppose different scatterers of the same source differ by one phase delay and a random amplitude gain. *α*_{il} is complex gain of the *il*th scatterer reflecting the reflection coefficient. *E*(*α*_{il}) = *α*_{i}. *n*_{mk}(*t*) are noises received.

(*θ*, *ϕ*; **u**_{i}) denotes the deterministic angular signal distribution function (ASDF) of the *i*th 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 **u**_{i} = (*θ*_{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 *s*_{i}(*t*) = *iLα*_{i} is the reflected signal of the *i*th 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 *i*th CD source can be written aswhere *s*_{1}(*t*) reflects of time behavior of 1st source and *ρ*_{i} is the complex correlation coefficient between the *i*th 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 *M*_{0} and *K*_{0} sensors along the direction of the *x* and *z* axes, respectively. Thus, subarrays smooth *M*_{s} = *M* − *M*_{0} + 1 times along the direction of the *x* axis and *K*_{s} = *K* − *K*_{0} + 1 times along the direction of the *z* axis.

Considering the (1, 1)th subarray, the receive vector of *M*_{0} 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 *k*th 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

From equations (10)–(12) we can conclude that the received vector of (*m*, *k*)th subarray along *x* axis can be obtained 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 *K*_{0} sensors along the *z* axis can be expressed aswhere is the generalized steering vector of , which can be defined as

According to equations (10)–(12), the received vector of closest *m*th and parallel to *z* axis array can be expressed aswhere **Φ**_{x} is the rotation operator, which can be written 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

From equations (10)–(12), the received vector of (*m*, *k*)th subarray along *x* axis can be obtained as follows:

Define as a vector selecting from 1 to *M*_{0}*K*_{0} − *M*_{0} row of . is a vector selecting from *M*_{0} + 1 to *M*_{0}*K*_{0} row of . and can be expressed aswhere (1:, *M*_{0}*K*_{0} − *M*_{0}:) denotes vectors from 1 to *M*_{0}*K*_{0} − *M*_{0} row of . is noise vector selecting from 1 to *M*_{0}*K*_{0} − *M*_{0} row of . is noise vector selecting from *M*_{0} + 1 to *M*_{0}*K*_{0} 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 **N**_{1}, **N**_{2}, **N**_{3,} and **N**_{4} 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 [31] has proved the spatial smoothing covariance matrices are full rank matrices when *M*_{0} > *q*, *K*_{0} > *q* and *M*_{S} ≥ *q*, *K*_{S} ≥ *q*.

Construct matrix **B** aswhere **B**_{1} is a *q* × *q* dimensional nonsingular matrix and **B**_{2} is a (4*M*_{0}*K*_{0} − 4*M*_{0} − *q*) × *q* matrix. There exists a *q* × (4*M*_{0}*K*_{0} − 4*M*_{0} − *q*) dimensional propagator operator **P** satisfying .

Letwhere **I**_{q×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 (4*M*_{0}*K*_{0} − 4*M*_{0}) × (*M*_{0}*K*_{0} − *K*_{0}) dimensional matrix into *q* × (*M*_{0}*K*_{0} − *K*_{0}) dimensional matrix **G** and (4*M*_{0}*K*_{0} − 4*M*_{0} − *q*) × (*M*_{0}*K*_{0} − *K*_{0}) matrix **H**:

The propagator operator **P** can be obtained aswhere (·)^{+} denotes the pseudoinverse. Divide **E** into four (*M*_{0}*K*_{0} − *M*_{0}) × *q* dimensional matrices **E**_{1}, **E**_{2}, **E**_{3,} and **E**_{4}.

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 *i*th 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 **E**_{1}, **E**_{2}, **E**_{3}, **E**_{4} 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*[(2*M*_{0}*K*_{0} − 2*M*_{0})^{3}], and eigendecomposition of **Ω**_{1} and **Ω**_{2}*O*(*q*^{3}). 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*[(2*MK* − 2*M*)^{3}], and eigendecomposition *O*(*q*^{3}). 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 *i*th 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 *i*th 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 *M*_{s} = *K*_{s} = 4, subarray elements parameter is *M*_{0} = *K*_{0} = 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 [18] which uses *L*-shaped arrays containing sensors in *x* and *z* axes, and the ESPRIT [23] 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°]. *M*_{s} = *K*_{s} = 4 and *M*_{0} = *K*_{0} = 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 [18] and ESPRIT [23] 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.

**(a)**

**(b)**

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. *M*_{s} = *K*_{s} = 4 and *M*_{0} = *K*_{0} = 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, *M*_{s} is supposed equal to *K*_{s,} and *M*_{0} is supposed equal to *K*_{0}. Figure 6 displays the performance of the proposed method versus the subarray number parameter *M*_{s} with *M*_{0} = 4, while Figure 7 shows the performance versus subarray elements parameter *M*_{0} with *M*_{s} = 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 *M*_{0} > *q*, *K*_{0} > *q*, *M*_{S} ≥ *q*, and *K*_{S} ≥ *q* and provides better performance as subarray number or subarray elements increases, but RMSE changes slightly when *M*_{0} or *M*_{s} 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 [30] 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. *M*_{s} = *K*_{s} = 3 and *M*_{0} = *K*_{0} = 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.

**(a)**

**(b)**

#### 5. Conclusions

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.

#### Appendix

#### A. Integral Expression of Received Signal

Considering equation (1), angle (*θ*_{il}, *ϕ*_{il}) which is DOA of *il*th scatterer of the *i*th 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 *i*th source are subject to the density distribution ,

Then, *y*_{mk}(*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

#### Data Availability

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.

#### Authors’ Contributions

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.

#### Acknowledgments

This research was funded by the National Natural Science Foundation of China (Grant numbers 61471299 and 51776222).