Optimization of array geometry for direction-of-arrival estimation using a priori information

This paper focuses on the estimation of the direction-of-arrival (DOA) of signals impinging on a sensor array. A novel method of array geometry optimization is presented that improves the DOA estimation performance compared to the standard uniform linear array (ULA) with half wavelength element spacing. Typically, array optimization only affects the beam pattern of a specific steering direction. In this work, the proposed objective function incorporates, on the one hand, a priori knowledge about the signal’s DOA in terms of a probability density function. By this means, the array can be adjusted to external conditions. On the other hand, a modified beam pattern expression that is valid for all possible signal directions is taken into account. By controlling the side lobe level and the beam width of this new function, DOA ambiguities, which lead to large DOA estimation errors, can be avoided. In addition, the DOA fine error variance is minimized. Using a globally convergent evolution strategy, the geometry optimization provides array geometries that significantly outperform the standard ULA with respect to DOA estimation performance. To show the quality of the algorithm, four optimum geometries are presented. Their DOA mean squared error is evaluated using the well known deterministic Maximum Likelihood estimator and compared to the standard ULA and theoretical lower bounds.


Introduction
This paper deals with the problem of estimating the direction-of-arrival (DOA) of signals impinging on an array of spatially distributed sensors.There exists a vast amount of DOA estimators, whose performance and accuracy compared to theoretical lower bounds like the Cramér-Rao lower Correspondence to: O. Lange (oliver.lange@de.bosch.com)bound have been intensively investigated.A good introduction to array signal processing and DOA estimation can be found in (Krim and Viberg, 1996) and the references therein.In this work, we focus on the deterministic Maximum Likelihood (DML) DOA estimator, which is asymptotically consistent and statistically efficient under certain regularity conditions (Stoica and Nehorai, 1989).Nevertheless, any other estimator that works with an arbitrary array geometry can be applied to the results of this work, as well.
Given a fixed number of sensor elements, the question of an optimum sensor placement with respect to DOA estimation performance naturally arises.In this context, it is intuitively clear that no "globally optimal" geometry for DOA estimation exists.Instead, an array geometry is always optimal with respect to certain presuppositions like the number of signals (Gershman and Böhme, 1997), the statistical properties of the DOA or the probability of large DOA estimation errors, etc.
Although there exist many different approaches to achieving an optimum geometry, basically all of them can be attributed directly or indirectly to the array beam pattern in terms of the side lobe level (SLL) and the shape of the main lobe.Because shape is a quantity that is difficult to characterize mathematically, the half power beam width (HPBW) is often used as a feature instead.It corresponds to the accuracy of the DOA estimates, i.e. the narrower the main beam, the lower the variance of the DOA estimates.Thus, the main beam width is directly related to the Cramér-Rao lower bound, which represents the minimum variance of an unbiased estimator (van Trees, 2002).Furthermore, a narrow beam width increases the possibility of angular signal separability in the multiple signal case.
It is well known in DOA estimation, that the mean squared error (MSE) departs from the Cramér-Rao lower bound when the signal-to-noise (SNR) ratio (or the sample size) falls below a specific limit.Beginning at this threshold, DOA estimates are dominated by large estimation errors, outliers whose probability of occurrence depends on the SLL of the beam pattern (Athley, 2005).
Moreover, in the multiple signal case, a low SLL prevents high SNR signals from overlapping low SNR signals.
By further reducing the SNR (or the sample size), the MSE goes to saturation, which can be approximated by the variance of a uniform distribution over the DOA search space (Richmond, 2006).Figure 1 quantitatively illustrates this typical DOA MSE curve and presents three regions of operation: 1.The asymptotic region is influenced by the shape of the main lobe, i.e. the HPBW; it can be approximated by the Cramër-Rao lower bound or the equivalent asymptotic variance of the estimator, respectively.
2. The threshold region is dominated by outliers, whose probability of occurrence is proportional to the SLL.
3. In the no information region, the SNR (or the sample size) is very low which leads to DOA estimates that are uniformly distributed over the search space.
All attempts of array geometry optimization that aim at DOA performance improvement have to take this MSE characteristic, which is based on a trade-off between SLL and beam width of the beam pattern, into account.In doing so, it should also be noted, that the beam pattern has to be unambiguous for all steering directions, where targets are to be expected.Ambiguities and high side lobes in this region should be avoided, as they lead to outliers in DOA estimation, even in high SNR regimes.As the usual beam pattern only represents the array response to a unit wave from a certain direction, we introduce a modified beam pattern (MBP), using a simple parameter transformation.By this means, the MBP allows for control of the HPBW and the SLL simultaneously for all valid steering directions.To take "real world" external conditions into account, we further incorporate a probability density function (PDF), defined for all valid steering directions.It is affected, for instance, by the angular a priori distribution of the target's DOA, the element factor of the sensors, the transmit beam pattern and the effect of blinds or radomes.This a priori PDF is applied as an angular dependent weighting factor to the MBP.By this means, we expand an approach by Oktel and Moses (2005), how previous knowledge about the DOA and hardware characteristics can easily be transferred into the array design process.
Using a globally convergent optimization routine, the HPBW and the SLL of the weighted MBP can now be controlled resulting in array geometries that outperform the standard uniform linear array (ULA) with half wavelength element spacing in terms of DOA estimation accuracy.
This paper is organized as follows.The data model and the MBP are introduced in Sect. 2. We also derive the weighted MBP, which additionally uses a priori knowledge in terms of a beta distributed PDF.Furthermore, we introduce the deterministic Maximum Likelihood DOA estimator and show, how prior knowledge can improve DOA estimation.In Sect.3, the objective function for array geometry optimization is defined.Some examples for optimum array geometries are presented and their DOA estimation performance is compared to the standard ULA using Monte Carlo simulations and the Cramér-Rao lower bound.

Data model and modified beam pattern
We consider a planar array geometry with N identical, omnidirectional sensors placed on a straight line, e.g. the x-axis (see Fig. 2).Q far-field sources that are orientated coplanar to the array are impinging on the array from azimuth directions 0 = [ 0,1 ,..., 0,Q ].The narrow band assumption is assumed to hold so that the complex baseband array output vector can be modeled by (van Trees, 2002) x where A(u 0 ) = a(u 0,1 ),...,a(u 0,Q ) is the (N × Q) array steering matrix with u 0,q = sin 0,q , q = 1,...,Q.The complex baseband source signals are denoted by s(k) = s 1 (k),...,s Q (k) T where (•) T means transpose, n(k) is the additive noise vector based on a Gaussian random process with zero mean and variance σ 2 n .The number of snapshots is denoted by K.The SNR of the q-th signal at each of the sensors is defined as The steering matrix A(u 0 ) consists of Q steering vectors a(u 0 ) whose n-th element is defined by Here, the element positions are denoted by p x,n and λ is the wavelength.

Derivation of the modified beam pattern
The normalized array beam pattern R(u,u 0 ) is defined as the squared magnitude of the response to a unit wave (s(k) = 1 ∀k) from direction u 0 in the case of no noise, i.e.
It can be clearly seen from Fig. 3, that knowledge of with −1 ≤ ũ ≤ 1, which represents the diagonal u 0 = −u in Fig. 3, allows for perfect reconstruction of all beam patterns 1 This means, that R(u,u 0 ) can be reduced to r( ũ) without any loss of information.Note that (Eq.6) is only valid for omnidirectional sensor elements.Non uniform element characteristics are taken into account in the next subsection where prior knowledge is considered.The function r( ũ) from (Eq. 5) is denoted as the modified beam pattern (MBP), because it can be obtained from the beam pattern by a simple rearrangement of the arguments u and u 0 .Figure 4 plots the MBP for the same parameters as in Fig. 3.Note that the region of unambiguous DOA estimates can be identified as

Definition of the a priori PDF
In many array signal processing applications, restrictive assumptions can be made concerning the statistical properties of the target's DOA.Hence, the probability of appearance or detection of a target at a certain direction u 0 is not uniform for all possible target directions −1 ≤ u ≤ 1. Often, a target is more likely to appear near the array boresight than laterally.
Or there might be a certain region, e.g. in surveillance applications, where the possibility of occurrence of a target is near zero.In addition, incoming signals from certain DOA's can be damped by the use of hardware constraints like blinds or radomes or by the influence of the array element factor.Furthermore, there are applications like automotive radar, where targets are actively illuminated by a transmit antenna.Neglecting the influence of multi path scattering, the DOA angular distribution directly corresponds to the beam shape of the transmit antenna.Thus, the mentioned characteristics directly or indirectly affect the statistical properties of the target's DOA.To account for the mentioned effects, we introduce the the PDF p u 0 (u) which represents the angular a priori information.In this work, we focus (without loss of generality) on the beta distribution where where a = b is chosen to assure a symmetric distribution with mean µ u 0 = 0 and variance σ 2 u 0 = 1 4(2a+1) .Figure 5 plots the PDF from (Eq. 8) for several values of u b and a.For a = 1, the PDF is uniform, i.e., there is no prior information about the DOA u 0 , except its restriction to the interval [−u b ,u b ].As a increases, the PDF becomes narrower.
To be consistent with the standard parameter estimation theory, we still consider the target's DOA u 0 as a deterministic parameter.Namely, we regard it as a realization of a stochastic process with the PDF from (Eq. 8).This is an important assumption if theoretical bounds like the Cramér-Rao lower bound on DOA estimation are considered.

Advantage of prior knowledge for DOA estimation: weighted deterministic Maximum Likelihood estimator
The derivation of the deterministic Maximum Likelihood (DML) DOA estimator can be found, for in (van Trees, 2002).In this work, we focus the single signal case, i.e. the estimation of u 0 .Hence, the DML estimate of u 0 is given by with the likelihood function where tr[•] is the trace operator and P ⊥ A is the projection matrix P ⊥ A = I − A(u) A H (u)A(u) −1 A H (u) onto the orthogonal complement of the column space of the steering matrix A(u).Furthermore, is the estimated spatial correlation matrix.
Any prior knowledge about the target's DOA statistical properties can improve the DML DOA estimation.By simweighting the likelihood function from (Eq. 10) with a priori PDF from (Eq. 8), we obtain the weighted istic Maximum Likelihood (WDML) estimate û0,WDML = arg max L −α (u)p u 0 (u) .
Here, the real valued parameter α ≥ 0 is used to control the ratio of the prior and the data: the larger the value of α, the lower the influence of the prior on the estimate.Note that (Eq.12) is the Maximum a posteriori estimator, if α represents the amount of collected data, i.e. the product of snapshots and channels α = KN (van Trees and Bell, 2007).
We will demonstrate this advantageous behavior with a quantitative example (see Fig. 6 for a graphical interpretation).Again, we refer to the 8-element ULA with element spacing 3 4 λ and assume one impinging source from u 0 = 1 2 .Due to the ambiguity at u = − 5 6 , the DML estimator (see dashed curve in Fig. 6) will probably provide a large estimation error.If, in addition, we assume, that the target's DOA follows a beta distribution as in (Eq.8) with a = 3, the WDML estimator (see solid curves in Fig. 6) reduces the effect of the ambiguity and provides the correct DOA estimate û0,WDML = 1 2 .

Optimization of the array geometry
In array geometry optimization, the objective function in general exhibits a multi modal character, i.e., local optimization routines might get stuck in local minima.Therefore, the optimization of the array geometry in this work uses the evolution strategy.A good survey of global optimization can be found in (Weise, 2007).Evolution strategy belongs to the family of evolutionary algorithms, which are based on biology-inspired methods like mutation, crossover, selection and survival of the fittest.In contrast to bit-encoded genetic algorithms, evolution strategies explore a real valued parameter space.The aim of this array geometry optimization is to identify a better array geometry concerning DOA estimation for one target, underlying a given a priori DOA PDF, compared to the standard ULA with half wavelength element spacing.

Definition of the objective function
In the array geometry optimization task, an objective function f (p x ) : R N → R has to be identified.For this purpose, we refer to the MBP (Eq.5) and the target's DOA PDF (Eq.8).We denote their product as weighted modified beam pattern (WMBP).To reduce the fine error variance, i.e. the Cramér-Rao lower bound, the main beam width has to be minimized.Therefore, the associated objective function f (p x ) calculates the HPBW of the WMBP from (Eq. 13): which is evaluated numerically.To avoid an increased possibility of gross errors, care is taken to keep the SLL of the WMBP below a threshold SLL thr .Thus, the array geometry optimization results in the constrained minimization problem min Figure 7 exemplarily illustrates the calculation of the objective function for the already mentioned 8-element ULA with sensor spacing 3 4 λ.At first, the product of the MBP and the PDF provides the WMBP, following (Eq.13).Here, we use the beta PDF from (Eq. 8) with u b = 0.5 and a = 1.Note that the SLL of the WMBP (SLL = −12.8dB) is below the threshold SLL thr = −10dB.Therefore, the HPBW of the WMBP is calculated using equation (Eq.14) with f (p x ) = ũ0,r − ũ0,l .

Optimization results
Three exemplary results of the geometry optimization will now be presented in order to illustrate the accuracy of the proposed objective function in conjunction with the evolution strategy.The DOA estimation performance for one signal in terms of the mean squared error (MSE) is evaluated for each generated geometry and compared to the standard ULA.The target's DOA u 0 is chosen to underlie the beta PDF, whose parameters are specified in the respective example.For DOA estimation, we use the WDML estimator from (Eq. 12) with α = 5, i.e., a high influence of the prior compared to the Maximum a posteriori estimator is considered.Each of the MSE curves is based on 6•10 4 Monte Carlo runs per SNR simulation point.The number of snapshots is K = 10 and we consider array geometries with N = 8 sensor elements.The respective optimized array geometries, including the standard ULA, are shown in Fig. 8.
In addition, the MSE performance plots also include the asymptotic variance  of the DML estimation error with the signal-to-noise ratio defined in (Eq.2) and the variance of the element positions (17) (Athley, 2005) and(van Trees, 2002).Note that in the single signal case, (Eq.16) is identical to the Cramér-Rao lower bound for a stochastic signal s(k) (Stoica and Nehorai, 1990).

Example 1
We choose the beta PDF with u b = 1 and a = 1, i.e., the DOA estimation should be unambiguous for the half space −1 ≤ u ≤ 1.In addition, the SLL threshold is set to SLL thr = −10dB.The optimized geometry is shown in Fig. 8b). Figure 9 plots the WMBP of the standard ULA and the optimized array.It can be seen that the ULA provides (near) ambiguous DOA estimates for | ũ| > 0.9.The optimized array, however, is unambiguous for the complete half space and exhibits a reduced HPBW at the expense of a slightly increased SLL.The associated DOA MSE curve is presented in the bottom of Fig. 9.It can be clearly seen, that in the case of u b = 1, i.e. a uniform distributed DOA u 0 ∈ [−1,1], the optimized array significantly outperforms the ULA due to the effect of outliers.If the ULA's ambiguous regions are not taken into account, i.e. u b = 0.9, its asymptotic MSE improves, but still does not reach the optimized array due to the marginal difference of the beam widths.Furthermore, Fig. 9 also shows, that the WDML estimator asymptotically approaches the asymptotic variance of the DML estimator, as it has been defined in (Eq.16).

Example 2
In this example, we assume that the signal's DOA is limited to −0.5 ≤ u 0 ≤ 0.5.Therefore, we choose the beta PDF with u b = 0.5 and a = 1 (see Fig. 5).To avoid ambiguities inside this region, the SLL is bounded above to SLL thr = −12dB and SLL thr = −8dB, respectively.The optimized geometries are shown in Fig. 8c) and d). Figure 10 plots the WMBP's for the ULA and the two optimized arrays.Due to the reduced DOA region, ambiguities outside of −0.5 ≤ u 0 ≤ 0.5 can now be accepted, as they have no impact on DOA estimation.This allows for an increased array aperture which leads to a smaller HPBW.The MSE curves in Fig. 10 show an improvement in terms of the SNR in the asymptotic region of 5dB and 10dB, respectively, compared to the ULA.

Example 3
Here, we assume a non-flat PDF with u b = 0.75 and a = 1.2, which has already been plotted in Fig. 5. Thus, in contrast to the other examples, the distribution of the DOA u 0 ∈ [−u b ,u b ] is no longer uniform.Again, the optimized array (see Fig. 8e) achieves a smaller HPBW at the expense of an increased SLL (see Fig. 11).Thus, the asymptotic MSE of the optimized array is reduced by 5.5dB.

Conclusion
In this paper, a novel method of array geometry optimization with the aim of DOA estimation improvement has been presented.Based on a modification of the standard beam pattern and an angular a priori density function, an objective function has been derived.By controlling the beam width and the SLL of the modified beam pattern, the MSE curve of DOA estimation is affected, while ambiguities are avoided.By the use of the a priori PDF, we have demonstrated, how effects from the "real world" like the DOA's statistical properties, antenna characteristics like the element factor and the transmit beam pattern or hardware constraints like radomes or blinds can easily be taken into account in the optimization process.Using an evolution strategy, have presented four optimum array geometries and evaluated their DOA estimation performance.Compared to the standard ULA, it has been shown, that the MSE can be significantly improved.

Fig. 1 .
Fig.1.Typical MSE curve of DOA estimation with three different regions of operation.In addition, likelihood functions are plotted corresponding to the asymptotic, the threshold and the no information region, respectively.

Fig. 2 .
Fig. 2.Array geometry with N elements on the x-axis with nonuniform inter element spacings and one far field source at azimuth angle .
Figure 3 exemplarily shows the beam patterns corresponding to

Fig. 5 .
Fig. 5. Different beta probability densities, transformed to u space, from (Eq. 8) with b = a for several values of u b and a.
1 dv is the beta function.The parameters a and b are real positive constants.A transformation of the random variable

Fig. 9 .
Fig.9.Top: WMBP for example 1 for the standard ULA and the optimized array.Bottom: DOA MSE performance of the WDML estimator for example 1 for a single signal.The DOA u 0 follows the PDF from (Eq. 8) with u b = 1 and u b = 0.9, respectively, and a = 1.

Fig. 10 .
Fig.10.Top: WMBP for example 2 for the standard ULA and the optimized arrays.Bottom: DOA MSE performance of the WDML estimator for example 2 for a single signal.The DOA u 0 follows the PDF from (Eq. 8) with u b = 0.5 and a = 1.

Fig. 11 .
Fig. 11.Top: WMBP's for example 3 for the standard ULA and the optimized array.Bottom: DOA MSE performance of the WDML estimator for example 3 for a single signal.The DOA u 0 follows the PDF from (Eq. 8) with u b = 0.75 and a = 1.2.