Adrien Bartoli

Piecewise Planar Segmentation for Automatic Scene Modeling

In this paper, we investigate the problem of the automatic creation of 3D models of man-made environments that we represent as collections of textured planes.

A typical approach is to automatically compute a sparse feature reconstruction and to manually give their planememberships as well as the delineation of the planes. Textures are then extracted from the images while optimizing the model, typically the disparity between marked and predicted edges.

We propose a means to automatically estimate the model of the scene, in terms of the number of planes and their parameters from a point feature reconstruction. The method is based on random sampling of reconstructed points to generate plane hypotheses. Each of these is then evaluated using a measure of approximate photoconsistency while recovering the corresponding plane delineation.

We then compute the maximum likelihood estimate of all scene parameters, i.e. the set of planes and reconstructed points as well as relative camera pose, with respect to actual images. The approach is validated on simulated data and real images.

Introduction

Automatic modeling of 3D scenes from images is one of the most challenging research areas in computer vision. Among the variety of all possible scenes, we consider the widespread case of man-made environments. Those scenes are often composed of piecewise planar primitives and can be modeled by a collection of planes [START_REF] Baillard | Automatic reconstruction of piecewise planar models from multiple views[END_REF][START_REF] Bartoli | Constrained structure and motion from N views of a piecewise planar scene[END_REF][START_REF] Debevec | Modeling and rendering architecture from photographs: a hybrid geometry-and image-based approach[END_REF][START_REF] Dick | Combining single view recognition and multiple view stereo for architectural scenes[END_REF] to a reasonable degree of approximation. Figure 1a shows an image of such a scene, overlaid by a piecewise planar segmentation, while figure 1b shows a view of the corresponding recovered model.

Existing systems are most of the time semi-automatic and can be seen as three-stage processes [START_REF] Debevec | Modeling and rendering architecture from photographs: a hybrid geometry-and image-based approach[END_REF][START_REF] Lang | 3D-city modeling with a digital one-eye stereo system[END_REF][START_REF] Streilein | Integration of digital photogrammetry and CAAD: Constraint-based modeling and semi-automatic measurement[END_REF]. First, a sparse 3D reconstruction of features (points, lines, etc.) as well as cameras is performed automatically [START_REF] Beardsley | 3D model acquisition from extended image sequences[END_REF]. From This work was supported by the project IST-1999-10756, VISIRE. there, two stages remain: choosing the scene model and estimating its parameters. The first stage is achieved by clustering reconstructed features into higher level geometric primitives such as cubes by e.g. marking edges in actual images. The second one consists in optimizing the quality of the model parameters, e.g. in minimizing the disparity between marked and predicted edges.

This approach has proven to give highly photorealistic results, but becomes effort-prone as the scene considered grows in complexity. In this paper, we devise a method that (a) (b) is automatic and produces the most likely scene parameters with respect to actual images.

In [START_REF] Morris | Image-consistent surface triangulation[END_REF], the scene surface is modeled as a set of triangles. The most likely triangulation with respect to actual images is computed using edge swaps from an initial solution obtained based on a Delaunay triangulation. The process is not guaranteed to converge to the global optimum and its complexity is high. The piecewise planarity is not taken into account, which reduces photorealism when rendering planar surfaces from oblique points of view.

To overcome these problems, we represent a scene as a collection of planes. This reduces the complexity of the model computation as well as its rendering and yields more photorealistic view synthesis of planar surfaces. The recovered structure and motion can be greatly enhanced by using the coplanarity information [START_REF] Bartoli | Constrained structure and motion from N views of a piecewise planar scene[END_REF]. We follow the same three stages described above. Our modifications are two-fold. First, we fill the gap of interactivity by automatically choosing a piecewise planar model of the scene and second, we compute the maximum likelihood estimate of the model parameters with respect to actual images.

Choosing the scene model can be done by fitting planes to reconstructed points. This is difficult because coplanar configurations that do not correspond to any world surface often arise in practice and have to be disambiguated using a physically meaningful criterion, i.e. photometric information. For that purpose, we propose to use a random sampling technique to hypothesize multiple plane equations. We then select the most likely planes with respect to actual images while checking for them to be distinct. Such an algorithm allows a point to lie on several planes.

We then refine the model parameters by non-linearly minimizing the squared greylevel difference between the actual and predicted images. This is the maximum likelihood estimate with respect to actual images. During this stage, we enforce the multi-coplanarity constraints. This is therefore a constrained photometric bundle adjustment.

This paper is organized as follow. In § §2 and 3, we review existing work and give our scene model and notations. In §4, we formally state the maximum likelihood estimation problem. We then propose an algorithm to choose the model of the scene in §5. Finally, we validate the method in §6, using both simulated data and real images.

Existing Work

Existing work on piecewise planar segmentation can be divided into two subsets, whether the photometric information, i.e. dense greylevel data, is taken into account or not. In particular, purely geometric criteria are sometimes used, see e.g. [START_REF] Alon | An integrated approach for segmentation and estimation of planar structure[END_REF][START_REF] Sinclair | Quantitative planar region detection[END_REF]. Such criteria do not allow to eliminate coplanar configurations that do not correspond to physical planes.

Photometric information is taken into account in [START_REF] Baillard | Automatic reconstruction of piecewise planar models from multiple views[END_REF][START_REF] Dick | Combining single view recognition and multiple view stereo for architectural scenes[END_REF]. The method described in [START_REF] Baillard | Automatic reconstruction of piecewise planar models from multiple views[END_REF] follows an approach based on line matching to segment planes from aerial images of urban scenes. In [START_REF] Dick | Combining single view recognition and multiple view stereo for architectural scenes[END_REF], the segmentation is performed using a collection of architectural primitives (pillars, doors, windows, etc.) with prior probabilities depending on the scene type.

Other possibilities for 3D modeling are space-carving [START_REF] Kutulakos | Approximate N-view stereo[END_REF] and image-consistent triangulation [START_REF] Morris | Image-consistent surface triangulation[END_REF], as described in the introduction.

Notations and Scene Model

We denote a scene as S = S m ∪ S p . As said before, it is essential to make a clear distinction between the scene model S m and its parameters S p .

The scene model S m = {l, θ m 1 , . . . , θ m l } includes the number of planes l and the plane models while the scene parameters S p = {θ p 1 , . . . , θ p l , Q, P} comprise the plane parameters, the cloud of 3D reconstructed points Q and the projection operator P, equivalent to reconstructed cameras. Similarly, each scene plane π is modeled by θ = θ m ∪θ p , i.e. it has a model and an associated set of parameters, both described in table 1. In more detail, V designates the set of images where π is visible and Π, Π and ∂ are defined respectively as two sets of points lying on the plane with respect to two different criteria described in §5 and the polygonal delineation of the plane, expressed as a list of points.

The set of n actual images is denoted as I. The projection operator P allows, given a scene S, to predict its images Î.

Notations #Π and conv Π designate the number of elements and the convexe envelope of a set of points Π.

Maximum Likelihood Estimation

In this section, we describe how to estimate the most likely parameters S p of the scene with respect to actual images, given the scene model S m .

If we assume our prior probability on the scene to be uniform, this is achieved by finding [START_REF] Dick | Combining single view recognition and multiple view stereo for architectural scenes[END_REF][START_REF] Morris | Image-consistent surface triangulation[END_REF]:

argmax S p

Pr(I|S),

where Pr(I|S) denotes the likelihood of the scene parameters S p for the scene model S m and the image set I. We derive an expression for this likelihood, in terms of the actual and predicted images.

We assume each pixel luminance to be corrupted by an i.i.d. centered gaussian noise. Using the same reasoning as in [START_REF] Morris | Image-consistent surface triangulation[END_REF], we obtain that the most likely parameters correspond to the minimum squared difference between the actual and predicted images.

The problem is then to find argmin S p C(S), where the cost function is defined by the RMS (Root Mean of Squares) error:

C(S) = ||I -Î|| 2 , (1)
where Î are the predicted images. This can be decomposed over the set of planes as:

C(S) = π∈S C(S π ),
where S π designates the scene restricted to plane π.

Given the scene model, the maximum likelihood estimate of its parameters can then be obtained by non-linearly minimizing (1) while enforcing the multi-coplanarity constraints between points and planes.

For that purpose, we minimally parameterize the estimation process using the technique given in [START_REF] Bartoli | Constrained structure and motion from N views of a piecewise planar scene[END_REF] so that multicoplanarity constraints are exactly enforced during the optimization. It is important to note that the cost function is not the same as that given in [START_REF] Bartoli | Constrained structure and motion from N views of a piecewise planar scene[END_REF], which is purely geometric, but similar to that given in [START_REF] Dick | Combining single view recognition and multiple view stereo for architectural scenes[END_REF]. The optimization is conducted using the Levenberg-Marquardt algorithm [START_REF] Gill | Practical Optimization[END_REF].

The process can be speeded up by estimating the error around the interest points only, as described in [START_REF] Baillard | Automatic reconstruction of piecewise planar models from multiple views[END_REF].

Piecewise Planar Segmentation

In this section, we present a method to compute a piecewise planar model of the scene, i.e. to choose S m .

Overview of the Algorithm

The algorithm inputs are the actual images as well as a cloud of reconstructed 3D points and camera matrices or equivalently the projection operator, see e.g. [START_REF] Beardsley | 3D model acquisition from extended image sequences[END_REF].

The piecewise planar segmentation is done by iteratively selecting the most likely plane using a random sampling technique. We use a multiple hypotheses version of RANSAC [START_REF] Fischler | Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography[END_REF] modified in two ways. First, as described in [START_REF] Torr | Robust computation and parametrization of multiple view relation[END_REF], we maximize the likelihood of the plane instead of its support (the support of a plane is the set of points lying on it, up to a certain tolerance). Secondly, we devise a segmentation scheme inspired from [START_REF] Alon | An integrated approach for segmentation and estimation of planar structure[END_REF] that allows for overlapping data segmentation, which is important in the piecewise planar segmentation of real scenes since a great number of planes are defined by points lying on several other planes.

Overlapping Data Segmentation

In this section, we present the robust estimator RANSAC and its application to disjoint data segmentation. We then propose a scheme that handles overlapping clusters.

RANSAC. This robust hypothesize and verify algorithm proceeds by repeatedly generating solutions estimated from minimal sets of points gathered from the data. This is the random sampling phase. It then tests each solution for support from the complete set of data. The parameters that maximize the support, i.e. the set of inliers, is then selected.

In the case of plane estimation, three points are randomly selected and a plane equation π is computed while testing for degeneracy. The set of points Π that geometrically lie on π, up to a predefined tolerance, is then computed imagebased [START_REF] Alon | An integrated approach for segmentation and estimation of planar structure[END_REF][START_REF] Sinclair | Quantitative planar region detection[END_REF] and the operation repeated. The dominant plane is the one that maximizes #Π. The number of iterations is computed as indicated in e.g. [START_REF] Fischler | Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography[END_REF].

As proposed in [START_REF] Torr | Robust computation and parametrization of multiple view relation[END_REF] in the case of image points, it is possible to use another cost function instead of #Π such as the likelihood Pr(S π |I) of the plane hypothesis. Computing this likelihood from the plane equation is described in §5. [START_REF] Bartoli | Constrained structure and motion from N views of a piecewise planar scene[END_REF].

Once the dominant plane has been estimated, points can be classified whether they lie, or not on the plane. In the former case, they belong to Π and in the latter, they are outliers.

Disjoint Data Segmentation. In the case of data segmentation, outliers have to be interpreted as points that do not lie on the dominant plane. Consequently, it is possible to recursively apply RANSAC on the set of outliers to perform data segmentation. This way of segmenting has proven to be efficient in the case of motion segmentation [START_REF] Demirdjian | Motion-egomotion discrimination and motion segmentation from image-pair streams[END_REF][START_REF] Torr | Outlier detection and motion segmentation[END_REF] or particular planar grouping [START_REF] Schaffalitzky | Planar grouping for automatic detection of vanishing lines and points[END_REF]. Indeed, in these cases the clusters are disjoint, e.g. a point will not, in the general case, satisfy two estimated different rigid motions.

In the case of piecewise planar segmentation, we can not make such an assumption of disjoint clusters since a point can lie on several planes. In the following paragraph, we present another scheme for data segmentation that allows overlapping clusters.

Overlapping Data Segmentation. In this case, we can not recursively apply RANSAC on the set of outliers for the reason mentioned above. Instead, we modify it to formulate multiple hypotheses, as described in [START_REF] Alon | An integrated approach for segmentation and estimation of planar structure[END_REF], so that a point can be selected as an inlier for more than one plane. The problem is that the same plane may be hypothesized more than once. This can be solved by constraining the random sampling phase to select a plane that has not been detected so far. Therefore, we have to define a measure between two planes π i and π j given their support, Π i and Π j . This measure must be physically meaningful, e.g. comparing the equations of the two planes will not give satisfactory results since it is a purely algebraic measure. We compare the sets of inliers Π i and Π j of the two planes. The two planes are identical if D(i, j) < γ where D(i, j) = #(Πi∩Πj ) #Πi+#Πj , i.e. the ratio of common inliers, and γ is a predefined rate, chosen as 0.4 < γ < 0.7 in our experiments, see §6. Another solution would be to assess each support with the other plane equation.

We have shown how to formulate plane equation hypotheses with respect to possibly overlapping clusters. In the following section, we examine how to estimate the likelihood of an hypothesis.

Retrieving Plane Parameters

So far, we have an hypothesis of a plane equation π as well as its geometric support Π. In order to evaluate the likelihood of this hypothesis, we have to determine its model θ m and parameters θ p (see table 1). that (Π ⊂ Π) ⇒ (∂ ⊂ conv Π). therefore, we compute conv Π and iteratively cut it to reach ∂.

Visibility images

For that purpose, we first warp each image of the plane with delineation conv Π, so that image points are aligned, and perform a Delaunay triangulation. The warping is done using inter-image plane homographies and bicubic interpolation. If we assume that the triangulation does not cut the edges of ∂, removing triangles from conv Π leads to ∂. Such a condition can be ensured by using an imageconsistent triangulation [START_REF] Morris | Image-consistent surface triangulation[END_REF]. Each triangle is then registered using an r-consistency check, as described in [START_REF] Kutulakos | Approximate N-view stereo[END_REF]. The dispersion radius r, expressed in pixels, controls the maximum image error of a reprojected point, allowing for non-perfectly planar surfaces and noisy images to be handled. The r-consistency check determines whether aligned triangles arising from different images represent the projection of the same planar surface, up to a tolerance modeled by r. The check fails if the greylevels of all pixels of all triangles are not shared by the other triangles in a disk of radius r around the pixel considered.

Such an algorithm allows for non convex and non-simply connected plane delineation recovery, provided an appropriate data structure in the latter case.

The set of points that lie inside ∂ is the global photometric support Π of the plane. Texture T . The texture is evaluated by warping each image of the plane onto a virtual fronto-parallel plane [START_REF] Liebowitz | Metric rectification for perspective images of planes[END_REF] and by taking the average of these values.

Experimental Results

In this section, we present the results we obtained by applying the algorithm described in this paper on both simulated data and real images.

Simulated Data

We compare our method to existing ones, notably to those consisting in using a disjoint data segmentation [START_REF] Dick | Combining single view recognition and multiple view stereo for architectural scenes[END_REF] or a purely geometric criterion [START_REF] Alon | An integrated approach for segmentation and estimation of planar structure[END_REF][START_REF] Sinclair | Quantitative planar region detection[END_REF].

The test bench consists of a cube of one meter side length observed by a set of cameras. Points are generated on the cube, possibly offset from their planes in order to simulate non-perfect coplanarity, and projected onto the images. Each face of the cube is texture mapped using random values. The texture is also projected onto the images and normalized to lie between 0 and 1. Up to 30, 5 and 1 points are generated on respectively each face, edge and vertex of the cube. Two cameras with a focal length of 1000 (expressed in number of pixels) and a 1 meter baseline are situated at a distance of 10 meters from the cube, such that 3 of its faces are observed. The intrinsic parameters are not supposed to be known which yields projective reconstructions.

In the sequel, we vary independently each of these parameters to compare the different approaches under various conditions. We measure the quality of segmentations using the absolute difference between the number of recovered planes and the number of simulated planes, i.e. 3. We use the median value over 100 trials. Matched image points are supposed to be known up to a gaussian centered noise of 1 pixel. The cloud of 3D points as well as camera matrices are obtained by minimizing reprojection errors in a bundle adjustment manner, using the ground truth as a starting point. The estimators compared are:

• DDS-geom: uses disjoint data segmentation and a purely geometric criterion, i.e. #Π, see §5.2;

• DDS-photo: idem but uses the photometric likelihood Pr(S π |I) criterion given in this paper in §5.3;

• ODS-geom: uses overlapping data segmentation described in §5.2 and a purely geometric criterion;

• ODS-photo: uses both the overlapping data segmentation and the photometric likelihood criterion. Let us describe the different experimental situations when varying a scene parameter and the simulation results we have obtained:

• Added image noise, figure 3a: a gaussian centered noise with standard deviation between 0.0 and 0.02 (i.e. 2% of the maximum) is added to the greylevel of each pixel;

• Number of points, figure 3b: the number of points is varied from 0 to 40 for each face of the cube. We also project each vertex but no points are simulated on the edges;

• Plane unflatness, figure 3c: 3D points are offset from the plane they lie on by distances drawn from a normal distribution with standard deviation between 0.0 and 0.1 meters.

We observe that in the general case, methods -photo, based on photometric information perform better than methodsgeom, based on a purely geometric criterion. In particular, the method ODS-photo performs better than the others in all cases.

In the case of added image noise, we observe that methods -geom perform worse than methods -photo. This can be explained by the fact that purely geometric criteria can not differentiate real planes from coplanar configurations that do not correspond to any simulated plane and as the number of simulated points is relatively low, such configurations might occur.

When varying the number of points, we observe that for low numbers (only the vertices of the cubes are used at the begining of the graph), methods DDS-perform worse. This is due to the fact that planes can be detected only when taking into account multi-coplanarity. When the number of points increases, methods -geom perform worse for the same reason as that given in the case of image noise.

In the case of increasing plane unflatness, methodsgeom perform worse since the probability of coplanar configurations that do not correspond to any simulated plane increases. The performance of methods -photo is not dramatically affected.

Real Images

In this section, we present the modeling results we obtained using the images shown in figure 4. Similar results have been obtained with other images.

We insist on the difficulty to build a piecewise planar model from these images. Indeed, windows undergo a significant parallax form their facade and can not be completely modeled using the given points which lie only on the facades. Moreover, their appearance changes across viewpoints because of specular effects. Several planes are also very thin.

We describe the different steps followed to perform a complete reconstruction, from the images to the 3D textured model. The camera has been off-line calibrated using a calibration pattern. However, it is clear that this information is only exploited for visualizing the results and does not affect the behaviour of the scene modeling.

Structure and motion initialization. Image points are obtained by tracking by hand. We perform a partial recon- struction from two images and incrementally add the others to obtain the complete structure and motion. We then run a bundle adjustment to minimize the reprojection error.

Piecewise planarity estimation. We use the algorithm described in this paper to recover the piecewise planarity of the scene. We have tried with a large dispersion radius [START_REF] Kutulakos | Approximate N-view stereo[END_REF], see §5. that could have been modeled separately, figure 5c. In the case of a low dispersion radius, 0.2% of the image size, which roughly corresponds to 1 pixel, the recovered piecewise planarity is incomplete, due to the parallax of windows relatively to the wall, figure 5a. An in-between dispersion radius of 1% of the image size, i.e. about 5 pixels, gives a complete and photorealistic-looking piecewise planar model, figure 5b.

The results we obtained are therefore satisfactory. The method tolerates approximately planar surfaces and is able to detect and model small planes.

Maximum likelihood estimate. We use the technique described in [START_REF] Bartoli | Constrained structure and motion from N views of a piecewise planar scene[END_REF] to minimally parameterize the structure while enforcing previously recovered piecewise planarity. The optimization process is conducted as described in §4 using the scene model S m obtained with r=5 pixels (figure 5b). The sparse feature reconstruction is exactly piecewise planar after constrained optimization. Figure 6 shows texturemapped renderings of the recovered model from various points of view different from the original ones. The texture maps for windows have been extracted from one image each, to avoid the blur effect caused by the parallax. Another possibility to render such approximate planes is to use the View-Dependent Texture-Mapping described in [START_REF] Debevec | Modeling and rendering architecture from photographs: a hybrid geometry-and image-based approach[END_REF].

Quality assessment. We have performed several measures on the structure obtained before and after the constrained photometric bundle adjustment for the piecewise planarity corresponding to a 1 pixel dispersion radius (fig-

σ 1 σ 2 µ
point-based 0.0138 0.0419 0.0350 plane-based 0.0087 0.0302 0.0209 ure 5a). Two kinds of quantity are significant: length ratios and angles. Table 2 shows measures of such quantities. In this table, σ 1 and σ 2 are the variances of the length of respectively the 6 vertical edges and the 6 horizontal edges of equal length, whereas µ is the mean of 1 -2α i /π where α i are the measures of right angles. The values given in table 2 show that the metric reconstruction obtained with the method given in this paper is of superior quality than the one obtained with a traditionnal method based only on points. This is important since several planes of the scene have a uniform appearance and therefore do not greatly contribute to the photometric bundle adjustment.

Similar measures performed on the other reconstructions show that the lower the dispersion radius is, the better the reconstruction is. This is due to the fact that high dispersion radii tend to merge close planes (see e.g. the two left planes of figure 5a, merged in figures 5b and 5c). For that reason, the individual position of each plane is less accurate and the reconstruction quality lower. Consequently, there is a trade-off between the reconstruction quality and the surface approximation.

Conclusions, Perspectives

We have investigated the automatic modeling of scenes using planes. The method is based on segmenting a cloud of reconstructed 3D points into multi-coplanar groups. We use random sampling to generate multiple plane hypotheses and select the most likely planes with respect to actual images. The use of a photometric measure allows to distinguish coplanar point configurations that do not correspond to any world plane from real planes. The algorithm allows for overlapping data segmentation which permits to detect the majority of planes.

The maximum likelihood estimate of the model with respect to actual images is then computed using an appropriate structure parameterization that enforces multicoplanarity constraints.

The use of simulated data shows that our method performs better, in terms of the number of detected planes, compared to those based on a purely geometric criterion or using a disjoint data segmentation scheme. We have also validated the method using real images. This shows a good rendering quality and demonstrates how the difficulty of approximate planarity of real planes can be overcome using approximate photoconsistency, i.e. r-consistency.

However, there are still limitations and further research to pursue. Indeed, the method depends on a judicious selection of feature points on edges and corners. This difficulty can be tackled using a matching process relying on piecewise planarity. We also plan to incorporate a model of light variation in the photoconsistency measure to specifically handle specular surfaces such as windows that often arise in practice.

Figure 1 .

 1 Figure 1. Example of piecewise planar modeling: (a) shows an image overlaid with the automatically recovered piecewise planarity and (b) the model rendered from a different point of view.

Figure 2 .

 2 Figure 2. Evaluating the delineation ∂ of a plane.

Figure 3 .

 3 Figure 3. Comparison of the error in the number of detected planes for different approaches versus different scene parameters.

Figure 4 .

 4 Figure 4. 3 out of the 6 actual images of an architectural scene overlaid with features. Note the significant parallax of windows relatively to the wall.

Figure 5 .

 5 Figure 5. Recovered piecewise planarity and rendering of the corresponding model using successively 0.2, 1 and 2% of the image size (roughly 1, 5 and 10 pixels) for the dispersion radius r. In the first case, (a), planes are accurately modeled, see e.g. the two left planes. The holes correspond to windows that can not be entirely modeled and are not approximated by a plane when using such a low tolerance. For the second case, (b), the tolerance is sufficiently high, so that windows can be approximated by a plane. Note that the two left planes of (a) are now merged. In the last case, (c), two planes suffice to completely approximate the scene surface, which yields several artefacts on the rendering.

Figure 6 .

 6 Figure 6. Texture-mapped model obtained.

Table 1 . Modeling a plane π is achieved by a set θ = θ m ∪ θ p where θ m is the plane model and θ p its parameters.

 1

		V	visibility images
	model θ m	Π individual geometric support Π global photometric support
		∂	polygonal delineation
	parameters θ p π T	plane equation texture map

Table 2 . Metric measures on the initial euclidean re- construction (point-based) and on that obtained after the constrained photometric bundle adjustment (plane- based). The lower

 2

σ 1 , σ 2 and µ (see text) are, the better the reconstruction is.


Thang Dinh

Alexandre Vu
email: vu@lss.supelec.fr

Rémy Renaux
email: renaux@lss.supelec.fr

Sylvie Boyer

Marcos
email: marcos@lss.supelec.fr

Some results on the Weiss-Weinstein bound for conditional and unconditional signal models in array processing

Keywords:

In this paper, the Weiss-Weinstein bound is analyzed in the context of sources localization with a planar array of sensors. Both conditional and unconditional source signal models are studied. First, some results are given in the multiple sources context without specifying the structure of the steering matrix and of the noise covariance matrix. Moreover, the case of an uniform or Gaussian prior are analyzed. Second, these results are applied to the particular case of a single source for two kinds of array geometries: a non-uniform linear array (elevation only) and an arbitrary planar (azimuth and elevation) array.

Introduction

Sources localization problem has been widely investigated in the literature with many applications such as radar, sonar, medical imaging, etc. One of the objective is to estimate the direction-of-arrival (DOA) of the sources using an array of sensors.

In array processing, lower bounds on the mean square error are usually used as a benchmark to evaluate the ultimate performance of an estimator. There exist several lower bounds in the literature. Depending on the assumptions about the parameters of interest, there are three main kinds of lower bounds. When the parameters are assumed to be deterministic (unknown), the main lower bounds on the (local) mean square error used are the well known Cramér-Rao bound and the Barankin bound (more particularly their approximations [START_REF] Mcaulay | A useful form of the Barankin lower bound and its application to PPM threshold analysis[END_REF][2][3] [START_REF] Todros | General classes of performance lower bounds for parameter estimation -part I: non-Bayesian bounds for unbiased estimators[END_REF]). When the parameters are assumed to be random with a known prior distribution, these lower bounds on the global mean square error are called Bayesian bounds [START_REF] Van Trees | Bayesian Bounds for Parameter Estimation and Nonlinear Filtering/Tracking[END_REF]. Some typical families of Bayesian bounds are the Ziv-Zakai family [START_REF] Ziv | Some lower bounds on signal parameter estimation[END_REF][7] [START_REF] Bell | Extended Ziv-Zakaï lower bound for vector parameter estimation[END_REF] and the Weiss-Weinstein family [START_REF] Weiss | A lower bound on the mean square error in random parameter estimation[END_REF][10] [START_REF] Renaux | A fresh look at the Bayesian bounds of the Weiss-Weinstein family[END_REF] [START_REF] Todros | General classes of performance lower bounds for parameter estimation -part II: Bayesian bounds[END_REF]. Finally, when the parameter vector is made from both deterministic and random parameters, the so-called hybrid bounds have been developed [START_REF] Rockah | Array shape calibration using sources in unknown locations-part I: Far-field sources[END_REF][14] [START_REF] Bay | On the hybrid Cramér-Rao bound and its application to dynamical phase estimation[END_REF].

Since the DOA estimation is a non-linear problem, the outliers effect can appear and the estimators mean square error exhibits three distinct behaviors depending on the number of snapshots and/or on the signal to noise ratio(SNR) [START_REF] Van Trees | Detection, Estimation and Modulation Theory[END_REF]. At high SNR and/or for a high number of snapshots, i.e., in the asymptotic region, the outliers effect can be neglected and the ultimate performance are described by the (classical/Bayesian/hybrid) Cramér-Rao bound. However, when the SNR and/or the number of snapshots decrease, the outliers effect lead to a quick increase of the mean square error: this is the so-called threshold effect. In this region, the behavior of the lower bounds are not the same. Some bounds, generally called global bounds (Barankin, Ziv-Zakai, Weiss-Weinstein) can predict the threshold while the others, called local bounds, like the Cramér-Rao bound or the Bhattacharyya bound cannot. Finally, at low SNR and/or at low number of snapshots, i.e., in the no-information region, the deterministic bounds exceed the estimator mean square error due to the fact that they do not take into account the parameter support. On the contrary, the Bayesian bounds exploit the parameter prior information leading to a "real" lower bound on the global mean square error.

In this paper 1 , we are interested in the Weiss-Weinstein bounds which is known to be one of the tightest Bayesian bound with the bounds of the Ziv-Zakai family. We will study the two main source models used in the literature [START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF]: the unconditional (or stochastic) model where the source signals are assumed to be Gaussian and the conditional (or deterministic) model where the source signals are assumed to be deterministic. Surprisingly, in the context of array processing, while closed-form expressions of the Ziv-Zakai bound (more precisizely its extension by Bell et. al. [START_REF] Bell | Explicit Ziv-Zakaï lower bound for bearing estimation[END_REF]) were proposed around 15 years ago for the unconditional model, the results concerning the Weiss-Weinstein bound are, most of the time, only conducted by way of computations. Concerning the unconditional model, in [START_REF] Nohara | Application of the Weiss-Weinstein bound to a two dimensional antenna array[END_REF], the Weiss-Weinstein bound has been evaluated by way of computations and has been compared to the mean square error of the MUSIC algorithm and classical Beamforming using a particular 8 × 8 element array antenna. In [START_REF] Nguyen | Comparison of performance bounds for DOA estimation[END_REF], the authors have introduced a numerical comparison between the Bayesian Cramér-Rao bound, the Ziv-Zakai bound and the Weiss-Weinstein bound for DOA estimation. In [21], numerical computations of the Weiss-Weinstein bound to optimize sensor positions for non-uniform linear arrays have been presented. Again in the unconditional model context, in [START_REF] Xu | Bayesian bounds for matched-field parameter estimation[END_REF], by considering the matched-field estimation problem, the authors have derived a semi closed-form expression of a simplified version of the Weiss-Weinstein bound for the DOA estimation. Indeed, the integration over the prior probability density function was not performed.

The conditional model (with known waveforms) is studied only in [START_REF] Renaux | Weiss-Weinstein bound for data aided carrier estimation[END_REF], where a closed-form expression of the WWB is given in the simple case of spectral analysis and in [START_REF] Vu | Closed-form expression of the Weiss-Weinstein bound for 3D source localization: the conditional case[END_REF] which is a simplified version of the bound.

While the primary goal of this paper is to give closed-form expressions of the Weiss-Weinstein bound for the DOA estimation of a single source with an arbitrary planar array of sensors, under both conditional and unconditional source signal models, we also provide partial closed-form expressions of the bound which could be useful for other problems. First, we study the general Gaussian observation model with parameterized mean or parameterized covariance matrix. Indeed, one of the success of the Cramér-Rao is that, for this observation model, a closed-form expression of the Fisher information matrix is available: this is the socalled Slepian-Bang formula [START_REF] Kay | Fundamentals of Statistical Signal Processing: Estimation Theory[END_REF]. Such kind of formulas have been less investigated in the context of bounds tighter than the Cramér-Rao bound. Second, some results are given in the multiple sources context without specifying the structure of the steering matrix and of the noise covariance matrix. Finally, these results are applied to the particular case of a single source for two kinds of array geometries: the non-uniform linear array (elevation only) and the planar (azimuth and elevation) array. Consequently, the aim of this paper is also to provide a textbook of formulas which could be applied in other fields. The Weiss-Weinstein bound is known to depend on parameters called test points and other parameters generally denoted s i . One particularity of this paper in comparison with the previous works on the Weiss-Weinstein bound is that we do not use the assumption

s i = 1/2, ∀i.
This paper is organized as follows. Section 2 is devoted to the array processing observation model which will be used in the paper. In Section 3, a short background on the Weiss-Weinstein bound is presented and two general closed-form expressions which will be the cornerstone for our array processing problems are derived. In Section 4 we apply these general results to the array processing problem without specifying the structure of the steering matrix. In Section 5, we study the particular case of the non-uniform linear array and of the planar array for which we provide both closed-form expressions of the bound in the context of a single stationnary source in the far field area. Some simulation results are proposed in Section 6. Finally, Section 7 gives our conclusions.

Problem setup

In this section, the general observation model generally used in array signal processing is presented as well as the first different assumptions used in the remain of the paper. Particularly, the so-called conditional and unconditional source models are emphasized.

Observations model

We consider the classical scenario of an array with M sensors which receives N complex bandpass

signals s (t) = [s 1 (t) s 2 (t) • • • s N (t)]
T . The output of the array is a M × 1 complex vector y (t) which can be modelled as follows (see, e.g., [START_REF] Van Trees | Detection, Estimation and Modulation theory: Optimum Array Processing[END_REF] or [START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF])

y (t) = A (θ) s (t) + n (t) , t = 1, . . . , T, (1)
where T is the number of snapshots, where θ

= [θ 1 θ 2 • • • θ q ]
T is an unknown parameter vector of interest2 , where A (θ) is the so-called M × N steering matrix of the array response to the sources, and where the M × 1 random vector n (t) is an additive noise.

Assumptions

• The unknown parameters of interest are assumed to be random with an a priori probability density function p (θ i ) , i = 1, . . . , q. These random parameters are assumed to be statistically independent such that the a priori joint probability density function is p (θ) = q ∏ i=1 p (θ i ). Note that this assumption will be only used in Subsections 4.2 and 4.3. We also assume that the parameter space, denoted Θ, is a connected subset of R q (see [START_REF] Haim | A comment on the Weiss-Weinstein bound for constrained parameter sets[END_REF]).

• The noise vector is assumed to be complex Gaussian, statistically independent of the parameters, i.i.d., circular, with zero mean and known covariance matrix

E [ n (t) n H (t) ] = R n .
This assumption will be made more restrictive in Section 5 where it will be assumed that R n = σ 2 n I. In any case, R n is assumed to be a full rank matrix.

• The steering matrix A (θ) is assumed such that the observation model is identifiable. From Section 3 to Section 4, the structure of A (θ) is not specified in order to obtain the more general results.

• Concerning the source signals, two kinds of models have been investigated in the literature (see, e.g., [START_REF] Stoica | Performances study of conditional and unconditional direction of arrival estimation[END_REF] or [START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF]) and will be alternatively used in this paper.

-M 1 : Unconditional or stochastic model : s(t) is assumed to be a complex circular random vector, i.i.d., statistically independent of the noise, Gaussian with zero-mean and known covariance matrix E [ s (t) s H (t) ] = R s . Note that concerning the previous results on the Cramér-Rao bound available in the literature [START_REF] Stoica | Performances study of conditional and unconditional direction of arrival estimation[END_REF], the covariance matrix R s is assumed to be unknown. In this paper, we have made the simpler assumption that the covariance matrix R s is known. These assumptions have already been used for the calculation of bounds more complex than the Cramér-Rao bound (see, e.g., [START_REF] Xu | Bayesian bounds for matched-field parameter estimation[END_REF], [START_REF] Bell | Explicit Ziv-Zakaï lower bounds for bearing estimation using planar arrays[END_REF], [START_REF] Reuven | The use of the Barankin bound for determining the threshold SNR in estimating the bearing of a source in the presence of another[END_REF]).

-M 2 : Conditional or deterministic model : ∀t, s(t) is assumed to be deterministic known. Note that, under the conditional model assumption, the signal waveforms can be assumed either unknown or known. While the conditional observation model with unknown waveforms seems more challenging, the conditional model with known waveforms signals which will be used in this paper can be found in several applications such as in mobile telecommunication and radar (see e.g.

[31], [START_REF] Cedervall | Efficient maximum likelihood DOA estimation for signals with known waveforms in presence of multipath[END_REF], and [START_REF] Li | Computationally efficient angle estimation for signals with known waveforms[END_REF]).

Likelihood of the observations

Let R y = E [ (y (t) -E [y (t)]) (y (t) -E [y (t)]) H ]
be the covariance matrix of the observation vector y (t) . According to the aforementioned assumptions, it is easy to see that under M 1 , the observations y (t) are distributed as a complex circular Gaussian random vector with zero mean and covariance matrix

R y (θ) = A(θ)R s A H (θ) + R n while under M 2 ,
p (Y; θ) = 1 π M T |R y (θ)| T exp ( - T ∑ t=1 y (t) H R -1 y (θ)y (t) ) , (2)
where R y (θ) = A(θ)R s A H (θ) + R n and the likelihood under M 2 is given by

p (Y; θ) = 1 π M T |R n | T exp ( - T ∑ t=1 (y (t) -A (θ) s (t)) H R -1 n (y (t) -A (θ) s (t))
) .

(3)

Weiss-Weinstein bound: Generalities

In this Section, we first remind to the reader the structure of the Weiss-Weinstein bound on the mean square error and the assumptions used to compute this bound. Second, a general result about the Gaussian observation model with parameterized mean or parameterized covariance matrix, which, to the best of our knowledge, does not appear in the literature is presented. This result will be useful to study both the unconditional model M 1 and the conditional model M 2 in the next Section.

Background

The Weiss-Weinstein bound for a q × 1 real parameter vector θ is a q × q matrix denoted WWB and is given as follows [START_REF] Weinstein | A general class of lower bounds in parameter estimation[END_REF]

WWB = HG -1 H T , (4)
where the q × q matrix H = [h 1 h 2 . . . h q ] contains the so-called test-points h i , i = 1, . . . , q such that θ + h i ∈ Θ ∀h i . The k, l-element of the q × q matrix G is given by

{G} k,l = E [( L s k (Y; θ + h k , θ) -L 1-s k (Y; θ -h k , θ) ) ( L s l (Y; θ + h l , θ) -L 1-s l (Y; θ -h l , θ) )] E [L s k (Y; θ + h k , θ)] E [L s l (Y; θ + h l , θ)] , (5)
where the expectations are taken over the joint probability density function p (Y, θ) and where the function

L (Y; θ + h i , θ) is defined by L (Y; θ + h i , θ) = p(Y,θ+h i ) p(Y,θ) . The notation L s k (Y; θ + h k , θ) means that s k is the power of L (Y; θ + h i , θ). The elements s i are such that s i ∈ [0, 1], i = 1, . . . , q.
Note that we have the following order relation [START_REF] Weinstein | A general class of lower bounds in parameter estimation[END_REF] Cov

( θ) = E [ ( θ -θ ) ( θ -θ ) T ] ≽ WWB, (6)
where A ≽ B means that the matrix A -B is a semi-positive definite matrix and where Cov

( θ)
is the global (the expectation is taken over the joint pdf p (Y, θ)) mean square error of any estimator θ of the parameter vector θ. Finally, in order to obtain a tight bound, one has to maximize WWB over the testpoints h i and s i i = 1, . . . , q. Note that this maximization can be done by using the trace of HG -1 H T or with respect to the Loewner partial ordering [START_REF] La Rosa | Barankin-type lower bound on multiple change-point estimation[END_REF]. In this paper we will use the trace of HG -1 H T which is enough to obtain tight results.

A general result on the Weiss-Weinstein bound and its application to the Gaussian observation models

An analytical result on the Weiss-Weinstein bound which will be useful in the following derivations and which could be useful for other problems is derived in this part. Note that this result is independent of the parameter vector size q and of the considered observation model.

Let us denote Ω the observation space. By rewriting the elements of matrix G (see Eqn. ( 5)) involved in the Weiss-Weinstein bound, one obtains for the numerator denoted by

N {G} k,l , N {G} k,l = E [( L s k (Y; θ + h k , θ) -L 1-s k (Y; θ -h k , θ) ) ( L s l (Y; θ + h l , θ) -L 1-s l (Y; θ -h l , θ) )] = ∫ Θ ∫ Ω p s k (Y, θ + h k ) p s l (Y, θ + h l ) p s k +s l -1 (Y, θ) dYdθ+ ∫ Θ ∫ Ω p 1-s k (Y, θ -h k ) p 1-s l (Y, θ -h l ) p 1-s k -s l (Y, θ) dYdθ - ∫ Θ ∫ Ω p s k (Y, θ + h k ) p 1-s l (Y, θ -h l ) p s k -s l (Y, θ) dYdθ- ∫ Θ ∫ Ω p 1-s k (Y, θ -h k ) p s l (Y, θ + h l ) p s l -s k (Y, θ) dYdθ, (7)
and for the denominator denoted by

D {G} k,l , D {G} k,l = E [L s k (Y; θ + h k , θ)] E [L s l (Y; θ + h l , θ)] = ∫ Θ ∫ Ω p s k (Y, θ + h k ) p s k -1 (Y, θ) dYdθ ∫ Θ ∫ Ω p s l (Y, θ + h l ) p s l -1 (Y, θ) dYdθ. (8)
Let us now define a function η (α, β, u, v) as

η (α, β, u, v) = ∫ Θ ∫ Ω p α (Y, θ + u) p β (Y, θ + v) p α+β-1 (Y, θ) dYdθ, (9)
where (α, β) ∈ [0, 1] 2 and where (u, v) are two q × 1 vectors such that θ + u ∈ Θ and θ + v ∈ Θ. The notation p α (Y, θ + u) means that α is the power of p (Y, θ + u) . By identification, it is easy to see that

{G} k,l = η (s k , s l , h k , h l ) + η (1 -s k , 1 -s l , -h k , -h l ) -η (s k , 1 -s l , h k , -h l ) -η (1 -s k , s l , -h k , h l ) η (s k , 0, h k , 0) η (0, s l , 0, h l ) . (10)
Note that we choose the arbitrary notation D {G} k,l = η (s k , 0, h k , 0) η (0, s l , 0, h l ) for the denominator.

The notation

D {G} k,l = η (s k , 1, h k , 0) η (1, s l , 0, h l ) or, even, D {G} k,l = η (s k , 0, h k , v) η (0, s l , u, h l ) will
lead to the same result.

With Eqn. [START_REF] Rapoport | WeissWeinstein lower bounds for markovian systems. part I: Theory[END_REF], it is clear that the knowledge of η (α, β, u, v) for a particular problem leads to the Weiss-Weinstein bound (without the maximization procedure over the test-points and over the parameters s i ). Surprisingly, this simple expression is given in [START_REF] Weinstein | A general class of lower bounds in parameter estimation[END_REF] only for s i = 1 2 , ∀i and not for the general case.

Let us now detail this function η (α, β, u, v). The function η (α, β, u, v) can be rewritten as

η (α, β, u, v) = ∫ Θ p α (θ + u) p β (θ + v) p α+β-1 (θ) ∫ Ω p α (Y; θ + u) p β (Y; θ + v) p α+β-1 (Y; θ) dYdθ = ∫ Θ ήθ (α, β, u, v) p α (θ + u) p β (θ + v) p α+β-1 (θ) dθ, (11)
where we define

ήθ (α, β, u, v, θ) = ∫ Ω p α (Y; θ + u) p β (Y; θ + v) p α+β-1 (Y; θ) dY. ( 12
)
Our aim is to give the most general result. Consequently, we will focus only on ήθ (α, β, u, v) since the a priori probability density function depends on the considered problem.

An important remark pointed out in [START_REF] Haim | A comment on the Weiss-Weinstein bound for constrained parameter sets[END_REF] is that the integration for the parameter space is with respect to the region {θ : p (θ) > 0} . However, since the functions being integrated are p (θ) , p (θ + u) , and p (θ + v) , then the actual region of integration (where all the functions are positive) is the intersection of three regions,

{θ : p (θ) > 0} ∩ {θ : p (θ + u) > 0} ∩ {θ : p (θ + v) > 0} .
Note that, in order to simplify the notation we only use Θ throughout this paper but this remark will be useful and explicitly specified in Section 4.2.

Gaussian observation model with parameterized covariance matrix

One calls (circular, i.i.d.) Gaussian observation model with parameterized covariance matrix, a model such that the observations y (t) ∼ CN (0, R y (θ)) where θ are the parameters of interest. Note that M 1 is a special case of this model since the parameters of interest appear only in the covariance matrix of the observations which has the following particular structure R y (θ

) = A(θ)R s A H (θ) + R n . The closed-form expression of ήθ (α, β, u, v) is given by: ήθ (α, β, u, v) = |R y (θ)| T (α+β-1) |R y (θ + u)| T α |R y (θ + v)| T β αR -1 y (θ + u) + βR -1 y (θ + v) -(α + β -1)R -1 y (θ) T . ( 13
)
The proof is given in Appendix .1. Note that, similar expressions are given in [START_REF] Bell | Explicit Ziv-Zakaï lower bound for bearing estimation[END_REF] (Eqn. (B.15)) and [START_REF] Van Trees | Detection, Estimation and Modulation Theory: Radar-Sonar Signal Processing and Gaussian Signals in Noise[END_REF] (p. 67, Eqn. ( 52)) for the particular case where α = s and β = 1s.

Gaussian observation model with parameterized mean

One calls (circular, i.i.d.) Gaussian observation model with parameterized mean, a model such that the observations y (t) ∼ CN (f (θ) , R y ) where θ are the parameters of interest. Note that M 2 is a special case of this model since the parameters of interest appear only in the mean of the observations which has the following particular structure f t (θ) = A(θ)s (t) (and

R y = R n ). The closed-form expression of ήθ (α, β, u, v) is given in this case by ln ήθ (α, β, u, v) = - T ∑ t=1 α (1 -α) f H t (θ + u) R -1 y f t (θ + u) +β (1 -β) f H t (θ + v) R -1 y f t (θ + v) + (1 -α -β) (α + β) f H t (θ) R -1 y f t (θ) -2 Re { αβf H t (θ + u) R -1 y f t (θ + v) +α (1 -α -β) f H t (θ + u) R -1 y f t (θ) + β (1 -α -β) f H t (θ + v) R -1 y f t (θ) } , (14)
or equivalently by

ln ήθ (α, β, u, v) = - T ∑ t=1 α (1 -α -β) R -1/2 y (f t (θ + u) -f t (θ)) 2 + αβ R -1/2 y (f t (θ + u) -f t (θ + v)) 2 +β (1 -α -β) R -1/2 y (f t (θ + v) -f t (θ)) 2 . ( 15
)
The details are given in Appendix .2.

General application to array processing

In the previous Section, it has been shown that the Weiss-Weinstein bound computation (or, at least, the matrix G computation) is reduced to the knowledge of the function η (α, β, u, v) given by Eqn. [START_REF] Weiss | A lower bound on the mean square error in random parameter estimation[END_REF]. As one can see in Eqn. [START_REF] Rapoport | WeissWeinstein lower bounds for markovian systems. part I: Theory[END_REF], the elements of the matrix G depend on η (α, β, u, v) for particular values of α, β, u, and v. Consequently, the goal of this Section is to detail these particular functions for our model given by Eqn. [START_REF] Mcaulay | A useful form of the Barankin lower bound and its application to PPM threshold analysis[END_REF]. Since Eqn. ( 9) can be decomposed into a deterministic part (in the sense where ήθ (α, β, u, v)

(see Eqn. ( 12
)
) only depends on the likelihood function) and a Bayesian part (when we have to integrate ήθ (α, β, u, v) over the a priori probability density function of the parameters), we will first focus on the particular functions ήθ (α, β, u, v) by using the results of the previous Section on the Gaussian observation model with parameterized mean or covariance matrix. Second, we will detail the passage from ήθ (α, β, u, v)

to η (α, β, u, v) in the particular case where p (θ i ) is a uniform probability density function ∀i. Another result will also be given in the case of a Gaussian prior.

Analysis of ήθ (α, β, u, v)

We will now detail the particular functions ήθ (α, β, u, v) involved in the different elements of {G} k,l , k, l ∈ {1, q} 2 for both models M 1 and M 2 .

Unconditional observation model M 1

Under the unconditional model M 1 , by using Eqn. [START_REF] Rockah | Array shape calibration using sources in unknown locations-part I: Far-field sources[END_REF], one obtains straightforwardly the functions

ήθ (α, β, u, v) involved in the elements {G} k,l = {G} l,k                                    ήθ (s k , s l , h k , h l ) = |Ry(θ)| T (s k +s l -1) |Ry(θ+h k )| T s k |Ry(θ+h l )| T s l |skR -1 y (θ+h k )+s l R -1 y (θ+h l )-(s k +s l -1)R -1 y (θ)| T , ήθ (1 -s k , 1 -s l , -h k , -h l ) = |Ry(θ)| T (1-s k -s l ) |Ry(θ-h k )| T (s k -1) |Ry(θ-h l )| T (s l -1) |(1-sk)R -1 y (θ-h k )+(1-s l )R -1 y (θ-h l )-(1-s k -s l )R -1 y (θ)| T , ήθ (s k , 1 -s l , h k , -h l ) = |Ry(θ)| T (s k -s l ) |Ry(θ-h l )| T (s l -1) |Ry(θ+h k )| T s k |skR -1 y (θ+h k )+(1-s l )R -1 y (θ-h l )-(s k -s l )R -1 y (θ)| T , ήθ (1 -s k , s l , -h k , h l ) = |Ry(θ)| T (s l -s k ) |Ry(θ-h k )| T (s k -1) |Ry(θ+h l )| T s l |(1-sk)R -1 y (θ-h k )+s l R -1 y (θ+h l )-(s l -s k )R -1 y (θ)| T , ήθ (s k , 0, h k , 0) = |Ry(θ)| T (s k -1) |Ry(θ+h k )| T s k |skR -1 y (θ+h k )-(s k -1)R -1 y (θ)| T , ήθ (0, s l , 0, h l ) = |Ry(θ)| T (s l -1) |Ry(θ+h l )| T s l |slR -1 y (θ+h l )-(s l -1)R -1 y (θ)| T . (16)
The diagonal elements of G are obtained by letting k = l in the above equations.

Conditional observation model M 2

Under the conditional model M 2 , by using Eqn. ( 15) with f t (θ) = A (θ) s (t) and R y = R n one obtains straightforwardly the functions ήθ (α, β, u, v) involved in the elements

{G} k,l = {G} l,k                                  ln ήθ (s k , s l , h k , h l ) = s k (s k + s l -1) ζ θ (h k , 0) + s l (s k + s l -1) ζ θ (h l , 0) -s k s l ζ θ (h k , h l ) , ln ήθ (1 -s k , 1 -s l , -h k , -h l ) = (s k -1) (s k + s l -1) ζ θ (-h k , 0) + (s l -1) (s k + s l -1) ζ θ (-h l , 0) -(1 -s k ) (1 -s l ) ζ θ (-h k , -h l ) , ln ήθ (s k , 1 -s l , h k , -h l ) = s k (s k -s l ) ζ θ (h k , 0) + (1 -s l ) (s k -s l ) ζ θ (-h l , 0) + s k (s l -1) ζ θ (h k , -h l ) , ln ήθ (1 -s k , s l , -h k , h l ) = (s k -1) (s k -s l ) ζ θ (-h k , 0) + s l (s l -s k ) ζ θ (h l , 0) + (s k -1) s l ζ θ (-h k , h l ) , ln ήθ (s k , 0, h k , 0) = s k (s k -1) ζ θ (h k , 0) , ln ήθ (0, s l , 0, h l ) = s l (s l -1) ζ θ (h l , 0) , (17)
where we define

ζ θ (µ, ρ) = T ∑ t=1 R -1/2 n (A (θ + µ) -A (θ + ρ)) s (t) 2 . ( 18
)
The diagonal elements of G are obtained by letting k = l in the above equations. Note that, since we are working on matrix G, all the previously proposed results are made whatever the number of test-points.

Analysis of η (α, β, u, v) with a uniform prior

Of course, the analysis of η (α, β, u, v) given by Eqn. [START_REF] Renaux | A fresh look at the Bayesian bounds of the Weiss-Weinstein family[END_REF] can only be conducted by specifying the a priori probability density functions of the parameters. Consequently, the results provided here are very specific. However, note that, in general, this aspect is less emphasized in the literature where most of the authors give results without specifying the prior probability density functions and compute the rest of the bound numerically (see e.g., [START_REF] Xu | Bayesian bounds for matched-field parameter estimation[END_REF][20] [START_REF] Bell | Performance bounds in parameter estimation with application to bearing estimation[END_REF]).

We assume that all the parameters θ i have a uniform prior distribution over the interval [a i , b i ] and are statistically independent. We will also assume one test-point per parameter, otherwise there is no possibility to obtain (pseudo) closed-form expressions. Consequently, the matrix H is such that

H = Diag ([h 1 h 2 • • • h q ]) , (19)
and the vector h i , i = 1, . . . , q, takes the value h i at the i th row and zero elsewhere. So, in this analysis, the vector u takes the value u i at the i th row and zero elsewhere and the vector v takes the value v j at the j th row and zero elsewhere (of course, we can have i = j). Under these assumptions, η (α, β, u, v) can be rewritten

3 for i ̸ = j η (α, β, u, v) = ∫ Θ ήθ (α, β, u, v) p α (θ i + u i ) p β (θ j + v j ) p β (θ i ) p α (θ j ) p α+β-1 (θ i ) p α+β-1 (θ j ) q ∏ k=1 k̸ =i,k̸ =j p (θ k ) dθ = 1 q ∏ k=1 (b k -a k ) ∫ Θ q-2 ∫ Θj ∫ Θi ήθ (α, β, u, v) dθ i dθ j d (θ/ {θ i , θ j }) , (20)
where

Θ i =    [a i , b i -u i ] if u i > 0, [a i -u i , b i ] if u i < 0, and Θ j =    [a j , b j -v j ] if v j > 0, [a j -v j , b j ] if v j < 0, . For i = j, one can have v = ±u, then one obtains η (α, β, u, v = ±u) = ∫ Θ ήθ (α, β, u, v) p α (θ i + u i ) p β (θ i ± u i ) p α+β-1 (θ i ) q ∏ k=1 k̸ =i p (θ k ) dθ = 1 q ∏ k=1 (b k -a k ) ∫ Θ q-1 ∫ Θi ήθ (α, β, u, v = ±u) dθ i d (θ/ {θ i }) . (21)
In the last equation

, if v = -u, then Θ i =    [a i + u i , b i -u i ] if u i > 0, [a i -u i , b i + u i ] if u i < 0, , while, if v = u, then Θ i =    [a i , b i -u i ] if u i > 0, [a i -u i , b i ] if u i < 0, .
Depending on the structure of ήθ (α, β, u, v) , η (α, β, u, v) has to be computed numerically or a closedform expression can be found.

Another particular case which appears sometimes is when the function ήθ (α, β, u, v) does not depend on θ (see, [START_REF] Renaux | Weiss-Weinstein bound for data aided carrier estimation[END_REF] [START_REF] Bell | Explicit Ziv-Zakaï lower bounds for bearing estimation using planar arrays[END_REF] and Section 5 of this paper). In this case, ήθ (α, β, u, v) is denoted ή (α, β, u, v) and one obtains from Eqn. ( 20)

[5][8][18][20][21][27]
η (α, β, u, v) = ή (α, β, u, v) q ∏ k=1 (b k -a k )     q ∏ k=1 k̸ =i,k̸ =j ∫ b k a k dθ k     ∫ Θi dθ i ∫ Θj dθ j = (b i -a i -|u i |) (b j -a j -|v j |) (b i -a i ) (b j -a j ) ή (α, β, u, v) , (22)
and from Eqn. ( 21)

η (α, β, u, v = u) = (b i -a i -|u i |) (b i -a i ) ή (α, β, u, v) , (23)
and

η (α, β, u, v = -u) = (b i -a i -2 |u i |) (b i -a i ) ή (α, β, u, v) . ( 24
)
4.3. Analysis of η (α, β, u, v) with a Gaussian prior Finally, one can mention that if the prior is now assumed to be Gaussian, i.e.,

θ i ∼ N ( µ i , σ 2 i )
∀i and ήθ (α, β, u, v) does not depend on θ one obtains after a straightforward calculation

η (α, β, u, v) = ή (α, β, u, v) ∫ R p α (θ i + u i ) p α-1 (θ i ) dθ i ∫ R p β (θ j + v j ) p β-1 (θ j ) dθ j = ή (α, β, u, v) exp ( - 1 2
( α (1 -α) u 2 i σ 2 i + β (1 -β) v 2 j σ 2 j )) , (25)
η (α, β, u, v = u) = ή (α, β, u, v) ∫ R p α+β (θ i + u i ) p α+β-1 (θ i ) dθ i = ή (α, β, u, v) exp ( - (α + β) (1 -α -β) u 2 i 2σ 2 i ) , (26)
and

η (α, β, u, v = -u) = ή (α, β, u, v) ∫ R p α (θ i + u i ) p β (θ i -u i ) p α+β-1 (θ i ) dθ i = ή (α, β, u, v) exp ( - ( α + β -α 2 -β 2 + 2αβ ) u 2 i 2σ 2 i ) . (27)

Specific applications to array processing: DOA estimation

We now consider the application of the Weiss-Weinstein bound in the particular context of source localization. Indeed, until now, the structure of the steering matrix A (θ) for a particular problem has not been used in the proposed (semi) closed-form expressions. Consequently, these previous results can be applied to a large class of estimation problems such as far-field and near-field sources localization, passive localization with polarized array of sensors, or radar (known waveforms).

Here, we want to focus on the direction-of-arrival estimation of a single source in the far-field area with narrow-band signal. In this case, the steering matrix A (θ) becomes a steering vector denoted a (θ) (except for one preliminary result concerning the conditional model which will be given whatever the number of sources in Section 5.1.2). The structure of this vector will be specified by the analysis of two kinds of array geometry: the non-uniform linear array from which only one angle-of-arrival can be estimated (θ becomes a scalar) and the arbitrary planar array from which both azimuth and elevation can be estimated (θ becomes a 2 × 1 vector). In any cases, the array always consists of M identical, omnidirectional sensors. Both model M 1 and M 2 will be considered and the noise will be assumed spatially uncorrelated: R n = σ 2 n I. Since we focus on the single source scenario, the variance of the source signal s (t) is denoted σ 2 s for the model M 1 . The general structure of the i th element of the steering vector is as follows

{a (θ)} i = exp ( j 2π λ r T i θ ) , i = 1, . . . , M (28)
where θ represents the parameter vector, where λ denotes the wavelength, and where r i denotes the coordinate of the i th sensor position with respect to a given referential. In the following, r i will be a scalar or a 2 × 1 vector depending on the context (linear array or planar array).

Preliminary results

Since our analysis is now reduced to the single source case, we give here some other closed-form expressions which will be useful when we will detail the specific linear and planar arrays.

Unconditional observation model M 1

In order to detail the set of functions ήθ given by Eqn. ( 16), one has to find closed-form expressions of the determinant |R y (θ + u)| and of determinants having the following structure:

m 1 R -1 y (θ 1 ) + m 2 R -1 y (θ 2 ) with m 1 + m 2 = 1 or m 1 R -1 y (θ 1 ) + m 2 R -1 y (θ 2 ) + m 3 R -1 y (θ 3 ) with m 1 + m 2 + m 3 = 1. Under M 1 , the observation covariance matrix is now given by R y (θ) = σ 2 s a(θ)a H (θ) + σ 2 n I M . (29)
Concerning the calculation of |R y (θ + u)|, it is easy to find

|R y (θ + u)| = σ 2M n ( 1 + σ 2 s σ 2 n ∥a(θ + u)∥ 2 ) . (30)
Moreover, after calculation detailed in Appendix .3, one obtains for the other determinants

m 1 R -1 y (θ 1 ) + m 2 R -1 y (θ 2 ) = 1 (σ 2 n ) M ( 1 -ϕ 1 m 1 ∥a(θ 1 )∥ 2 + m 2 ϕ 2 ∥a(θ 2 )∥ 2 -ϕ 1 m 1 ϕ 2 m 2 ( a H (θ 1 )a(θ 2 ) 2 -∥a(θ 1 )∥ 2 ∥a(θ 2 )∥ 2 )) (31)
and

m 1 R -1 y (θ 1 ) + m 2 R -1 y (θ 2 ) + m 3 R -1 y (θ 3 ) = 1 (σ 2 n ) M ( 1 - 3 ∑ k=1 m k ϕ k ∥a(θ k )∥ 2 - 1 2 3 ∑ k=1 3 ∑ k ′ =1 k ′ ̸ =k m k ϕ k m k ′ ϕ k ′ ( a H (θ k )a(θ k ′ ) 2 -∥a(θ k )∥ 2 ∥a(θ k ′ )∥ 2 ) - ( 3 ∏ k=1 m k ϕ k )    3 ∏ k=1 ∥a(θ k )∥ 2 - 1 2 3 ∑ k=1 3 ∑ k ′ =1 k ′ ̸ =k 3 ∑ k ′′ =1 k ′′ ̸ =k ′ ̸ =k a H (θ k )a(θ k ′ ) 2 ∥a(θ k ′′ )∥ 2 +a H (θ 3 )a(θ 2 )a H (θ 1 )a(θ 3 )a H (θ 2 )a(θ 1 ) + a H (θ 3 )a(θ 1 )a H (θ 1 )a(θ 2 )a H (θ 2 )a(θ 3 ) )) , (32)
where

ϕ k = σ 2 s σ 2 s ∥a(θ k )∥ 2 + σ 2 n , k = 1, 2, 3. (33)

Conditional observation model M 2

Note that the results proposed here are in the context of any number of sources. Under the conditional model, the set of functions ήθ given by Eqn. ( 17) is linked to the function ζ θ (µ, ρ) given by Eqn. [START_REF] Bell | Explicit Ziv-Zakaï lower bound for bearing estimation[END_REF]. In this analysis, the vector µ takes the value µ i at the i th row and zero elsewhere and the vector ρ takes the value ρ j at the j th row and zero elsewhere (of course, one can has i = j). In Appendix .4, the calculation of the following closed-form expressions for ζ θ (µ, ρ) are detailed.

• If (m -1) p + 1 ≤ i, j ≤ mp, where p denotes the number of parameters per source, then, we have

ζ θ (µ, ρ) = T ∑ t=1 ∥{s (t)} m ∥ 2 M ∑ i=1 M ∑ j=1 { R -1 n } i,j exp ( j 2π λ ( r T j -r T i ) θ m ) × ( exp ( -j 2π λ r T i µ m ) -exp ( -j 2π λ r T i ρ m
)) ( exp

( j 2π λ r T j µ m ) -exp ( j 2π λ r T j ρ m )) . ( 34
)
• Otherwise, if (m -1) p + 1 ≤ i ≤ mp and (n -1) p + 1 ≤ j ≤ np, then we have

ζ θ (µ, ρ) = T ∑ t=1 ∥{s (t)} m ∥ 2 M ∑ i=1 M ∑ j=1 { R -1 n } i,j exp ( j 2π λ ( r T j -r T i ) θ m ) exp ( -j 2π λ r T i µ m ) exp ( j 2π λ r T j µ m ) + T ∑ t=1 ∥{s (t)} n ∥ 2 M ∑ i=1 M ∑ j=1 { R -1 n } i,j exp ( j 2π λ ( r T j -r T i ) θ n ) exp ( -j 2π λ r T i ρ n ) exp ( j 2π λ r T j ρ n ) -2 Re ( T ∑ t=1 {s (t)} * m {s (t)} n × M ∑ i=1 M ∑ j=1 { R -1 n } i,j exp ( j 2π λ ( r T j θ n -r T i θ m ) ) exp ( -j 2π λ r T i µ m ) exp ( j 2π λ r T j ρ n )   . ( 35
)
In particular, if one assumes R n = σ 2 n I, then, several simplifications can be done:

• If (m -1) p + 1 ≤ i, j ≤ mp, then ζ θ (µ, ρ) = 1 σ 2 n M ∑ i=1 exp ( -j 2π λ r T i µ m ) -exp ( -j 2π λ r T i ρ m ) 2 T ∑ t=1 ∥{s (t)} m ∥ 2 , ( 36
)
where we note that the function ζ θ (µ, ρ) does not depend on the parameter θ.

• Otherwise, if (m -1) p + 1 ≤ i ≤ mp and (n -

1) p + 1 ≤ j ≤ np, then ζ θ (µ, ρ) = 1 σ 2 n M ∑ i=1 exp ( -j 2π λ r T i µ m ) 2 T ∑ t=1 ∥{s (t)} m ∥ 2 + 1 σ 2 n M ∑ i=1 exp ( -j 2π λ r T i ρ n ) 2 T ∑ t=1 ∥{s (t)} n ∥ 2 -2 Re ( 1 σ 2 n M ∑ i=1 exp ( j 2π λ r T i (θ n -θ m ) ) exp ( -j 2π λ r T i µ m ) exp ( j 2π λ r T i ρ n ) T ∑ t=1 {s (t)} * m {s (t)} n ) ( 37
)
It is clear that the proposed above formulas for both the unconditional and the conditional models can be applied to any kind of array geometry and whatever the number of sources. However, they generally depend on the parameter vector θ. This means that, in general, the calculation of the set of functions η will have to be performed numerically (except if one is able to find a closed-form expression of Eqn. ( 11)). In the following we present a kind of array geometry where, fortunately, the set of functions ήθ will not depend on θ leading to a straightforward calculation of the bound.

3D Source localization with a planar array

We first consider the problem of DOA estimation of a single narrow band source in the far field area by using an arbitrary planar array. In fact, we start by this general setting because the non-uniform linear array is clearly a particular case of this array. Without loss of generality, we assume that the sensors of this array lay on the xOy plan with Cartesian coordinates (see Fig. .1). Therefore, the vector r i contains the coordinate of the i th sensor position with respect to this referential, i.e.,

r i = [d xi d yi ] T , i = 1, . . . , M .
From [START_REF] Stoica | Performances study of conditional and unconditional direction of arrival estimation[END_REF], the steering vector is given by

a(θ) = [ exp ( j 2π λ (d x1 u + d y 1 v) ) . . . exp ( j 2π λ (d xM u + d yM v) )] T , (38)
where, as in [START_REF] Bell | Explicit Ziv-Zakaï lower bound for bearing estimation[END_REF], the parameter vector of interest is θ = [u v] T where

   u = sin ϕ cos φ, v = sin ϕ sin φ, (39)
and where ϕ and φ represent the elevation and azimuth angles of the source, respectively. The parameters space is such that u ∈ [-1, 1] and v ∈ [-1, 1]. Therefore, we assume that they both follow a uniform distribution over [-1, 1]. Note that from a physical point of view, it should be more tempting to choose a uniform prior for ϕ and φ. This will lead to a probability density functions for u and v not uniform. To the best of our knowledge, this assumption has only been used in the context of lower bounds in [START_REF] Nguyen | Comparison of performance bounds for DOA estimation[END_REF]. Unfortunately, such prior leads to an untractable expression of the bound (see Eqn. (21) of [START_REF] Nguyen | Comparison of performance bounds for DOA estimation[END_REF]). Consequently, other authors have generally not specified the prior leading to semi closed-form expressions of bounds (i.e. that it remains a numerical integration to perform over the parameters) [START_REF] Nguyen | Comparison of performance bounds for DOA estimation[END_REF][37] [START_REF] Xu | Bayesian bounds for matched-field parameter estimation[END_REF]. On the other hand, in order to obtain a closed-form expression, authors have generally used a simplified assumption, i.e. a uniform prior directly on u and v (see, for example, [21] [START_REF] Xu | A bound on mean-square estimation error with background parameter mismatch[END_REF]). In this paper, we have followed the same way by expecting a slight modification of performance with respect to a more physical model and in order to be able to get closed-form expressions of the bound.

We choose the matrix of test points such that

H = [h u h v ] =   h u 0 0 h v   . (40)
Then, we have:

θ+h u = [u + h u v] T and θ+h v = [u v + h v ]
T . Moreover, we now have two elements

s i ∈ [0, 1] , i = 1
, 2 for which we will prefer the notation s u and s v , respectively.

Unconditional observation model M 1

Under M 1 , let us set

U SN R = σ 4 s σ 2 n (M σ 2 s +σ 2 n ) .
The closed-form expressions of the elements of matrix

G =   {G} uu {G} uv {G} vu {G} vv 
 are given by (see Appendix .5 for the proof):

{G} uu =           ( 1 -|hu| 2 ) ( 1 + 2su(1 -2su)USNR ( M 2 - M ∑ k=1 exp ( -j 2π λ dx k hu ) 2 )) -T + ( 1 -|hu| 2 ) ( 1 + 2(1 -su)(2su -1)USNR ( M 2 - M ∑ k=1 exp ( -j 2π λ dx k hu ) 2 )) -T -2 (1 -|hu|) ( 1 + su(1 -su)USNR ( M 2 - M ∑ k=1 exp ( -j 4π λ dx k hu ) 2 )) -T           ( 1 -|hu| 2 ) 2 ( 1 + s u (1 -s u )U SN R ( M 2 - M ∑ k=1 exp ( -j 2π λ d x k h u ) 2 )) -2T , (41)
{G} vv =           ( 1 -|hv | 2 ) ( 1 + 2sv(1 -2sv)USNR ( M 2 - M ∑ k=1 exp ( -j 2π λ dy k hv ) 2 )) -T + ( 1 -|hv | 2 ) ( 1 + 2(1 -sv)(2sv -1)USNR ( M 2 - M ∑ k=1 exp ( -j 2π λ dy k hv ) 2 )) -T -2 (1 -|hv|) ( 1 + sv(1 -sv)USNR ( M 2 - M ∑ k=1 exp ( -j 4π λ dy k hv ) 2 )) -T           ( 1 -|hv| 2 ) 2 ( 1 + s v (1 -s v )U SN R ( M 2 - M ∑ k=1 exp ( -j 2π λ d y k h v ) 2 )) -2T , ( 42
)       1 -U SN R           susv ( M ∑ k=1 exp ( -j 2π λ (dx k hu -dy k hv) ) 2 -M 2 ) +su(1 -su -sv) ( M ∑ k=1 exp ( -j 2π λ dx k hu ) 2 -M 2 ) +sv(1 -su -sv) ( M ∑ k=1 exp ( -j 2π λ dy k hv ) 2 -M 2 )           -susv(1 -su -sv) U 2 SN R σ 2 n σ 2 s × ×              M ∑ k=1 exp ( j 2πdy k hv λ ) M ∑ k=1 exp ( -j 2πdx k hu λ ) M ∑ k=1 exp ( j 2π(dx k hu -dy k hv ) λ ) + M ∑ k=1 exp ( -j 2πdy k hv λ ) M ∑ k=1 exp ( j 2πdx k hu λ ) M ∑ k=1 exp ( -j 2π(dx k hu-dy k hv ) λ ) -M M ∑ k=1 exp ( -j 2π λ dy k hv ) 2 -M M ∑ k=1 exp ( -j 2π λ dx k hu ) 2 -M M ∑ k=1 exp ( -j 2π λ (dx k hu -dy k hv)       (1 -su)(1 -sv) ( M ∑ k=1 exp ( j 2π λ (dx k hu -dy k hv) ) 2 -M 2 ) +(1 -su)(su + sv -1) ( M ∑ k=1 exp ( j 2π λ dx k hu ) 2 -M 2 ) +(1 -sv)(su + sv -1) ( M ∑ k=1 exp ( j 2π λ dy k hv ) 2 -M 2 )           -(1 -su)(1 -sv)(su + sv -1) U 2 SN R σ 2 n σ 2 s × ×              M ∑ k=1 exp ( j 2πdy k hv λ ) M ∑ k=1 exp ( -j 2πdx k hu λ ) M ∑ k=1 exp ( j 2π(dx k hu-dy k hv ) λ ) + M ∑ k=1 exp ( -j 2πdy k hv λ ) M ∑ k=1 exp ( j 2πdx k hu λ ) M ∑ k=1 exp ( -j 2π(dx k hu -dy k hv ) λ ) -M M ∑ k=1 exp ( -j 2π λ dy k hv ) 2 -M M ∑ k=1 exp ( -j 2π λ dx k hu ) 2 -M M ∑ k=1 exp ( -j 2π λ (dx k hu -dy k hv) 1 -USNR           su(1 -sv) ( M ∑ k=1 exp ( -j 2π λ (dx k hu + dy k hv) ) 2 -M 2 ) +su(sv -su) ( M ∑ k=1 exp ( -j 2π λ dx k hu ) 2 -M 2 ) +(1 -sv)(sv -su) ( M ∑ k=1 exp ( j 2π λ dy k hv ) 2 -M 2 )           -su(1 -sv)(sv -su) U 2 SN R σ 2 n σ 2 s × ×              M ∑ k=1 exp ( j 2πdy k hv λ ) M ∑ k=1 exp ( j 2πdx k hu λ ) M ∑ k=1 exp ( -j 2π(dx k hu+dy k hv ) λ ) + M ∑ k=1 exp ( -j 2πdy k hv λ ) M ∑ k=1 exp ( -j 2πdx k hu λ ) M ∑ k=1 exp ( j 2π(dx k hu +dy k hv ) λ ) -M M ∑ k=1 exp ( j 2π λ dy k hv ) 2 -M M ∑ k=1 exp ( -j 2π λ dx k hu ) 2 -M M ∑ k=1 exp ( -j 2π λ (dx k hu + dy k hv)       sv(1 -su) ( M ∑ k=1 exp ( -j 2π λ (dx k hu + dy k hv) ) 2 -M 2 ) +sv(su -sv) ( M ∑ k=1 exp ( -j 2π λ dx k hu ) 2 -M 2 ) +(1 -su)(su -sv) ( M ∑ k=1 exp ( -j 2π λ dy k hv ) 2 -M 2 )           -sv(1 -su)(su -sv) U 2 SN R σ 2 n σ 2 s × ×              M ∑ k=1 exp ( j 2πdy k hv λ ) M ∑ k=1 exp ( j 2πdx k hu λ ) M ∑ k=1 exp ( -j 2π(dx k hu+dy k hv ) λ ) + M ∑ k=1 exp ( -j 2πdy k hv λ ) M ∑ k=1 exp ( -j 2πdx k hu λ ) M ∑ k=1 exp ( j 2π(dx k hu +dy k hv ) λ ) -M M ∑ k=1 exp ( -j 2π λ dy k hv ) 2 -M M ∑ k=1 exp ( -j 2π λ dx k hu ) 2 -M M ∑ k=1 exp ( -j 2π λ (dx k hu + dy k hv) cos ( 2π λ dy k hv ) ) -2susvCSNR ( M - M ∑ k=1 cos ( 2π λ (dx k hu -dy k hv) ) )      + exp      2(su -1)(su + sv -1)CSNR ( M - M ∑ k=1 cos ( 2π λ dx k hu
) )

+2(sv -1)(su + sv -1)CSNR

( M - M ∑ k=1 cos ( 2π λ dy k hv ) ) -2(1 -su)(1 -sv)CSNR ( M - M ∑ k=1 cos ( 2π λ (dx k hu -dy k hv) ) )      -exp      2su(su -sv)CSNR ( M - M ∑ k=1 cos ( 2π λ dx k hu ) ) +2(1 -sv)(su -sv)CSNR ( M - M ∑ k=1 cos ( 2π λ dy k hv ) ) +2su(sv -1)CSNR ( M - M ∑ k=1 cos ( 2π λ (dx k hu + dy k hv) ) )      -exp      2(su -1)(su -sv)CSNR ( M - M ∑ k=1 cos ( 2π λ dx k hu ) ) +2sv(sv -su)CSNR ( M - M ∑ k=1 cos ( 2π λ dy k hv
) )

+2(su -1)svCSNR ( M -

Unconditional observation model M 1

The closed-form expression of the unconditional Weiss-Weinstein bound, denoted U W W B, is given by

U W W B = h 2 θ ( 1 -|h θ | 2 ) 2 ( 1 + s(1 -s)U SN R ( M 2 - M ∑ k=1 exp ( -j 2π λ d x k h θ ) 2 )) -2T             ( 1 -|h θ | 2 )       ( 1 + 2s(1 -2s)U SN R ( M 2 - M ∑ k=1 exp ( -j 2π λ d x k h θ ) 2 )) -T + ( 1 + 2(1 -s)(2s -1)U SN R ( M 2 - M ∑ k=1 exp ( -j 2π λ d x k h θ ) 2 )) -T       -2 (1 -|h θ |) ( 1 + s(1 -s)U SN R ( M 2 - M ∑ k=1 exp ( -j 4π λ d x k h θ ) 2 )) -T             . (55)
In order to find one optimal value of s that maximizes HG -1 H T , ∀h θ we have considered the derivative of HG -1 H T w.r.t. s. The calculation (not reported here) is straightforward and it is easy to see that

∂HG -1 H T ∂s s= 1 2 = 0.
Consequently, the Weiss-Weinstein bound has just to be optimized over h θ and is simplified leading to

U W W B = sup h θ h 2 θ ( 1 -|h θ | 2 ) 2 ( 1 + USNR 4 ( M 2 - M ∑ k=1 exp ( -j 2π λ d x k h θ ) 2 )) -2T 2 ( 1 -|h θ | 2 ) -2 (1 -|h θ |) ( 1 + USNR 4 ( M 2 - M ∑ k=1 exp ( -j 4π λ d x k h θ ) 2 )) -T . (56)
In the classical case of a uniform linear array (i.e., d x k = d), this expression can be still simplified by noticing that

M ∑ k=1 exp ( -j 2π λ d x k h θ ) =M exp ( -j 2πd λ h θ
) .

Conditional observation model M 2

The closed-form expression of the conditional Weiss-Weinstein bound CW W B is given by

CW W B = h 2 θ ( 1 -|h θ | 2 ) 2 exp ( 4s(s -1)C SN R ( M - M ∑ k=1 cos ( 2π λ d x k h θ ) ))         ( 1 -|h θ | 2 )     exp ( 4s(2s -1)C SN R ( M - M ∑ k=1 cos ( 2π λ d x k h θ ) )) + exp ( 4(2s -1)(s -1)C SN R ( M - M ∑ k=1 cos ( 2π λ d x k h θ ) ))     -2 (1 -|h θ |) exp ( 2s(s -1)C SN R ( M - M ∑ k=1 cos ( 4π λ d x k h θ ) ))         . (57)
Again, it is easy to check that ∂HG -1 H T ∂s s= 1 2 = 0. Consequently, one optimal value of s that maximizes

HG -1 H T , ∀h θ is s = 1 2 .
The Weiss-Weinstein bound is then simplified as follows

CW W B = sup h θ h 2 θ ( 1 -|h θ | 2 ) 2 exp ( -C SN R ( M - M ∑ k=1 cos ( 2π λ d x k h θ ) )) 2 ( 1 -|h θ | 2 ) -2 (1 -|h θ |) exp ( -1 2 C SN R ( M - M ∑ k=1 cos ( 4π λ d x k h θ ) )) . (58)
In the classical case of a uniform linear array (i.e., d x k = d), this expression can be still simplified by noticing that

M ∑ k=1 cos ( 2π λ d x k h θ ) = M cos ( 2πd λ h θ
) .

Simulation results and analysis

As an illustration of the previously derived results, we first consider the scenario proposed in Fig. 5 To the best of our knowledge, their are no closed-form expressions of the Ziv-Zakai bound for the conditional model available in the literature. In this case, we consider 3D source localization using a Vshaped array. Indeed, it has been shown that this kind of array is able to outperform other classical planar arrays, more particularly the uniform circular array [START_REF] Gazzah | Cramér-Rao bounds for antenna array design[END_REF]. This array is made from two branches of uniform linear arrays with 6 sensors located on each branches and one sensor located at the origin. We denote ∆ the angle between these two branches. The sensors are equally spaced with a half-wavelength. The number of snapshots is T = 20. Fig. .3 shows the behavior of the Weiss-Weinstein bound with respect to the opening angle ∆. One can observe that when ∆ varies, the estimation performance concerning the estimation of parameter u varies slightly. On the contrary, the estimation performance concerning the estimation of parameter v is strongly dependent on ∆. When ∆ increases from 0 • to 90 • , the Weiss-Weinstein bound of v decreases, as well as the SNR threshold. Fig. .3 also shows that ∆ = 90 • is the optimal value, which is different with the optimal value ∆ = 53.13 • in [START_REF] Gazzah | Cramér-Rao bounds for antenna array design[END_REF] since the assumptions concerning the source signal are not the same.

Conclusion

In this paper, the Weiss-Weinstein bound on the mean square error has been studied in the array processing context. In order to analyze the unconditional and conditional signal source models, the structure of the bound has been detailed for both Gaussian observation models with parameterized mean or parameterized covariance matrix.

1 σ 2 n ( I M - σ 2 s a(θ)a H (θ) σ 2 s ∥a(θ)∥ 2 + σ 2 n ) .
Then,

3 ∑ k=1 m k R -1 y (θ k ) = 1 σ 2 n 3 ∑ k=1 m k ( I - σ 2 s a(θ k )a H (θ k ) σ 2 s ∥a(θ k )∥ 2 + σ 2 n ) . ( .10)
Since the rank of a(θ k )a H (θ k ) is equal to 1 and since θ 1 ̸ = θ 2 ̸ = θ 3 (except for h k = h l = 0), the above matrix has M -3 eigenvalues equal to 1

σ 2 n 3 ∑ k=1
m k and 3 eigenvalues corresponding to the eigenvectors made from the linear combination of a(θ 1 ), a(θ 2 ), and a(θ 3 ): a(θ 1 ) + pa(θ 2 ) + qa(θ 3 ). The determinant will then be the product of these M eigenvalues5 . Let us set

ϕ k = σ 2 s σ 2 s ∥a(θ k )∥ 2 + σ 2 n , k = 1, 2, 3. ( .11)
Then, the three aforementioned eigenvalues denoted λ must satisfy:

( 3 ∑ k=1 m k R -1 y (θ k )
) (a(θ 1 ) + pa(θ 2 ) + qa(θ 3 )) = λ (a(θ 1 ) + pa(θ 2 ) + qa(θ 3 )) .

(.12) By using Eqn. (.10) in the above equation and after a factorization with respect to a(θ 1 ), a(θ 2 ), and a(θ 3 ) one obtains ( xm 1 ϕ 1 ∥a(θ 1 )∥ ) a(θ 2 ) + ( -m 3 ϕ 3 a H (θ 3 )a(θ 1 )m 3 ϕ 3 pa H (θ 3 )a(θ 2 ) + q ( xm 3 ϕ 3 ∥a(θ 3 )∥

2

)) a(θ 3 ) = 0, (

where 6

x = 1σ 2 n λ.

(

Consequently, the coefficients of a(θ 1 ), a(θ 2 ), and a(θ 3 ) are equals to leading to a system of three equations with two unknown (p and q). Solving the two first equations to find 7 p and q, and applying the solution into the last equation, one obtains the following polynomial equation of x

x 3 -x 2 3 ∑ k=1 m k ϕ k ∥a(θ k )∥ 2 - x 2 3 ∑ k=1 3 ∑ k ′ =1 k ′ ̸ =k m k ϕ k m k ′ ϕ k ′ ( a H (θ k )a(θ k ′ ) 2 -∥a(θ k )∥ 2 ∥a(θ k ′ )∥ 2 )
-m 1 m 2 m 3 ϕ 1 ϕ 2 ϕ 3 ( ∥a(θ 1 )∥ 2 ∥a(θ 2 )∥ 2 ∥a(θ 3 )∥ 2a H (θ 2 )a(θ 3 ) 2 ∥a(θ 1 )∥

2

a H (θ 1 )a(θ 2 ) 2 ∥a(θ 3 )∥ 2a H (θ 3 )a(θ 1 ) 2 a H (θ 2 ) 2 + a H (θ 3 )a(θ 2 )a H (θ 1 )a(θ 3 )a H (θ 2 )a(θ 1 ) +a H (θ 3 )a(θ 1 )a H (θ 1 )a(θ 2 )a H (θ 2 )a(θ 3 ) ) = 0

Since we are only interested by the product of the three eigenvalues, we do not have to solve this polynomial in λ and only the opposite of the last term is required. This leads to Eqn. ( 31) with

of [ 18 ]

 18 , i.e., the DOA estimation under the unconditional model using an uniform circular array consisting of M = 16 sensors with a half-wavelength inter-sensors spacing. The numbers of snapshots is T = 100. Since the array is symmetric, the performance estimation concerning the parameters u and v are the same, this is why only the performance with respect to the parameters u is given in Fig..2. The Weiss-Weinstein bound is computed using Eqn. (45), (46) and (47). The Ziv-Zakai bound is computed using Eqn.[START_REF] Vu | Closed-form expression of the Weiss-Weinstein bound for 3D source localization: the conditional case[END_REF] in[START_REF] Bell | Explicit Ziv-Zakaï lower bound for bearing estimation[END_REF]. The empirical global mean square error (MSE) of the maximum a posteriori (MAP) estimator is obtained over 2000 Monte Carlo trials. As in the Fig. (1b) in [18], one observes that both the Weiss-Weinstein bound and the Ziv-Zakai bound are tight w.r.t. the MSE of the MAP and capture the SNR threshold. Note that, in the Fig. (1b) in [18], the Weiss-Weinstein bound was computed numerically only.

3 ∑m k = 1 . 2 ∑

 312 k=1Of course, the closed-form expression of m 1 R -1 y (θ 1 ) + m 2 R -1 y (θ 2 ) is obtained by letting m 3 = 0 and k=1 m k = 1 in Eqn.[START_REF] Cedervall | Efficient maximum likelihood DOA estimation for signals with known waveforms in presence of multipath[END_REF].Appendix .4. Closed-form expressions of ζ θ (µ, ρ)Remind that the function ζ θ (µ, ρ) is defined by Eqn.[START_REF] Bell | Explicit Ziv-Zakaï lower bound for bearing estimation[END_REF]. Let us define p as the number of parameters per sources (assumed to be constant for each sources). Then, without loss of generality, the full parameter vector θ can be decomposed as θ =[ θ T 1 . . . θ T N ] T where θ i = [θ i,1 . . . θ i,p ]T , i = 1, . . . , N with q = N p.Remind that µ = [0 . . . µ i . . . 0] T and ρ = [ 0 . . . ρ j . . . 0 ] T . It exists two distinct cases to study: when both index i and j are such that (m -1)p + 1 ≤ i ≤ mp, m = 1, . . . , N and (m -1)p + 1 ≤ j ≤ mp or when6 Note that, from Eqn. (16),

 1 + 1 +Figure . 1 :

 111 Figure .1: 3D source localization using a planar array antenna.

Figure . 2 :

 2 Figure .2: Ziv-Zakai bound, Weiss-Weinstein bound and empirical MSE of the MAP estimator: unconditional case.

2 Figure . 3 :

 23 Figure .3: Weiss-Weinstein bounds of the V-shaped array w.r.t. the opening angle ∆.

  the observations y (t) are distributed as a complex circular Gaussian random vector with mean A(θ)s (t) and covariance matrix R y = R n . Moreover, in both case the

	observations are i.i.d..
	Therefore, the likelihood, p (Y; θ) , of the full observations matrix Y = [y (1) y (2) . . . y (T )] under M 1
	is given by

  2pm 1 ϕ 1 a H (θ 1 )a(θ 2 )qm 1 ϕ 1 a H (θ 1 )a(θ 3 ) ) a(θ 1 ) + ( -m 2 ϕ 2 a H (θ 2 )a(θ 1 ) + p ( xm 2 ϕ 2 ∥a(θ 2 )∥ 2 ) qm 2 ϕ 2 a H (θ 2 )a(θ 3 )

  m 1 ϕ 1 ∥a(θ 1 )∥ 2 ) ( xm 2 ϕ 2 ∥a(θ 2 )∥ 2 ) m 1 ϕ 1 m 2 ϕ 2 a H (θ 1 )a(θ 2 )a H (θ 2 )a(θ 1 ) m 1 ϕ 1 a H (θ 1 ) ( m 2 ϕ 2 a(θ 2 )a H (θ 2 ) + ( xm 2 ϕ 2 ∥a(θ 2 )∥ 2

	7 p and q are given by	M ∑ k=1	exp	(	j	2πdy k hv λ
				p =  m 2 ϕ 2 a H (θ 2 ) m 1 ϕ 1 a H (θ 1 )	( m 1 ϕ 1 a(θ 1 )a H (θ 1 ) + ( m 2 ϕ 2 a(θ 2 )a H (θ 2 ) +	( x -m 1 ϕ 1 ∥a(θ 1 )∥ 2 ( x -m 2 ϕ 2 ∥a(θ 2 )∥ 2	) )	I I	) )	a(θ 3 ) a(θ 3 )	,	(.15)
				
	and	q =	(	   x ) 	I	)	a(θ 3 )	.	(.16)

Section 5.2.2 of this paper has been partially presented in[START_REF] Vu | Closed-form expression of the Weiss-Weinstein bound for 3D source localization: the conditional case[END_REF]

Note that one source can be described by several parameters. Consequently, q > N in general.

In this case, one has to have a particular attention to the integration domain as mentionned in Section 3.2. It will not be the case for the Gaussian prior since the support is R.

Note that we are only interested by the eigenvalues. Consequently, the linear combination of of a(θ 1 ), a(θ 2 ), and a(θ 3 ) can be written a(θ 1 ) + pa(θ 2 ) + qa(θ 3 ) instead of ra(θ 1 ) + pa(θ 2 ) + qa(θ 3 )

and, of course, {G} uv = {G} vu . Consequently, the unconditional Weiss-Weinstein bound is 2 × 2 matrix given by:

which has to be optimized over s u , s v , h u , and h v . Concerning the optimization over s u and s v , several other works in the literature have suggested to simply use s u = s v = 1/2. Most of the time, numerical simulations of this simplified bound compared with the bound obtained after optimization over s u and s v leads to the same results while their is no formal proof of this fact (see [START_REF] Van Trees | Bayesian Bounds for Parameter Estimation and Nonlinear Filtering/Tracking[END_REF] page 41 footnote [START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF]. Note that, thanks to the expressions obtained in the next Section concerning the linear array, we will be able to prove that s = 1/2 is a (maybe not unique) correct choice for any linear array. In the case of the planar array treated in this Section, we will only check this property by simulation.

In the particular case where s u = s v = 1/2 one obtains the following simplified expressions

and

Again, the Weiss-Weinstein bound is obtained by using the above expressions in Eqn. (44) and after an optimization over the test points. The optimization over the test points can be done over a search grid or by using the ambiguity diagram of the array in order to reduce significantly the computational cost (see [START_REF] Reuven | A Barankin-type lower bound on the estimation error of a hybrid parameter vector[END_REF], [START_REF] Xu | Bayesian bounds for matched-field parameter estimation[END_REF], [START_REF] Reuven | The use of the Barankin bound for determining the threshold SNR in estimating the bearing of a source in the presence of another[END_REF], [START_REF] Tabrikian | Barankin bounds for source localization in an uncertain ocean environment[END_REF]).

Conditional observation model M 2

Under M 2 , let us set

The closed-form expressions of the elements of matrix G are given by (see Appendix .6 for the proof):

and {G} uv = {G} vu . Consequently, the conditional Weiss-Weinstein bound is 2 × 2 matrix given by using the above equations in Eqn. (44). As for the unconditional case, if we set s u = s v = 1/2, one obtains the following simplified expressions

By using the above expressions in Eqn. (44) and after an optimization over the test points, one obtains the Weiss-Weinstein bound.

Source localization with a non-uniform linear array

We now briefly consider the DOA estimation of a single narrow band source in the far area by using a non-uniform linear array antenna. Without loss of generality, let us assume that the linear array antenna lays on the Ox axis of the coordinate system (see Fig.

We assume that the parameter θ follows a uniform distribution over [-1, 1]. As in Section 4.2 and since the parameter of interest is a scalar, matrix H of the test points becomes a scalar denoted h θ . In the same way, there is only one element s i ∈ [0, 1] which will be simply denoted s. The closed-form expressions given here are straightforwardly obtained from the aforementioned results on the planar array about the element denoted {G} uu . We will continue to use the previously introduced notations

and

Appendix .1. Closed-form expression of ήθ (α, β, u, v) under the Gaussian observation model with parameterized covariance

Since y (t) ∼ CN (0, R y (θ)), one has,

where

one has

Appendix .2. Closed-form expression of ήθ (α, β, u, under the Gaussian observation model with parameterized mean

Let us set x = y -(αf t (θ + u)

And ξ (t) can be rewritten as

4 For simplicity, the dependance on t of f and y is not emphasized.

22

where

Note that ξ (t) is independent of x. By defining X = [x (1) , x (2) , . . . , x (T )], the function ήθ (α, β, u, v) becomes ήθ (α, β, u, = 1

since

Note that this calculation is actually an extension of the result obtained in Appendix A of [START_REF] Xu | Bayesian bounds for matched-field parameter estimation[END_REF] in which

2 and m 3 = 0, but follows the same method. The inverse of R y can be deduced from the Woodbury formula

Appendix .4.1. The case where (m -1) p + 1 ≤ i, j ≤ mp

In this case, one has:

and consequently,

Due to Eqn. ( 28), one has

In particular, in the case where R n = σ 2 n I one obtains

Appendix .4.2. The case where (m -1) p + 1 ≤ i ≤ mp and where (n -1) p + 1 ≤ j ≤ np Without loss generality, we assume that n > m. Then,

and consequently,

Let us set

By using the structure of the steering matrix A, it leads to (30) and (32) into Eqn. ( 16), and by considering the following expressions The set of functions ήθ (α, β, u, v) is given by Eqn. [START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF]. So, it only remains the calculation of functions

Again, since the set of functions ζ θ (µ, ρ) does not depend on θ, the set of functions ήθ (α, β, u, v) is given by plugging the above equations into Eqn. [START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF] and does not depend on θ. Consequently, as in unconditional case, the set of functions η (α, β, u, v) is obtained by using the results of Section 4.2 whatever the considered prior on θ. In our case of a uniform prior, the results are straightforward and leads to Eqn. (48), (49) and (50).