Lionel Bombrun
email: lionel.bombrun@gipsa-lab.inpg.fr

Jean-Marie Beaulieu
email: jean-marie.beaulieu@ift.ulaval.ca

SEGMENTATION OF POLARIMETRIC SAR DATA BASED ON THE FISHER DISTRIBUTION FOR TEXTURE MODELING

Keywords: Polarimetric SAR images, Fisher distribution, Texture, KummerU, Classification, Segmentation

The Polarimetric Synthetic Aperture Radar (PolSAR) covariance matrix is generally modeled by a complex Wishart distribution. For textured scenes, the product model is used and the texture component is often modeled by a Gamma distribution. In this paper, authors propose to use the Fisher distribution for texture modeling. From a Fisher distributed texture component, we derive the distribution of the complex covariance matrix and we propose to implement the KummerU distribution in a hierarchical segmentation and a hierarchical clustering algorithm. Segmentation and classification results are shown on synthetic images and on ESAR L-band PolSAR data over the Oberpfaffenhofen test-site.

INTRODUCTION

SAR data are the result of a coherent imaging system that produces the speckle noise phenomenon. If the antenna can emit and receive vertically and horizontally polarized waves, the covariance matrix is used to characterize the interactions between the electromagnetic wave and the scatterers. For fully developed speckle, the covariance matrix follows the complex Wishart distribution [START_REF] Goodman | Statistical analysis based on a certain multivariate complex Gaussian distribution (an introduction)[END_REF]. It is assumed that land cover backscatter characteristics are homogeneous (uniform or not textured). For textured scenes, the "scalar product model" has been proposed [START_REF] Joughin | Polarimetric Density Functions for Multilook Polarimetric Signatures[END_REF]. The observed signal is the product of a positive scalar texture component µ (polarization independent) with the speckle component. The texture term is generally modeled by a Gamma distribution and the covariance matrix follows a K distribution.

In this paper, we propose to use a more general probability density function (pdf) for the texture parameter: the Fisher distribution. Then, with the scalar product model, we show that the pdf of the observed covariance matrix can be expressed in term of the KummerU confluent hypergeometric function of the second kind. Then, we propose to implement this new pdf in a hierarchical segmentation algorithm. Segmentation results are illustrated on a synthetic image. Finally, the KummerU pdf is implemented in a totally unsupervised hierarchical classification algorithm. Classification results are shown on ESAR L-band PolSAR data over the Oberpfaffenhofen test-site.

COVARIANCE MATRIX DISTRIBUTION

Pdf for homogeneous scenes

One-look PolSAR data are completely characterized by the scattering matrix S = S hh S hv S vh S vv . S hv is the scattering component for horizontally polarized transmitting antenna and vertically polarized receiving antenna. For the monostatic case, the reciprocity assumption holds. The cross-polarization terms S hv and S vh are equal. The radar backscattering signal is described by the target scattering vector x = S hh , S vv , √ 2S hv T where T is the transposition operator. For homogeneous (no textured) scene (denoted by subscript h) and fully developed speckle, the pdf of x follows a zero mean multivariate complex Gaussian distribution [START_REF] Goodman | Statistical analysis based on a certain multivariate complex Gaussian distribution (an introduction)[END_REF]. The L-look covariance matrix,

Z h = (1/L) L k=1 x k x H k , follows the complex Wishart distribution [1] p Z h (Z h |Σ h ) = L Lp |Z h | L-p exp -L tr Σ h -1 Z h π p(p-1) 2 Γ(L) • • • Γ(L -p + 1)|Σ h | L (1)
where

Σ h = E xx H is the population covariance ma- trix. E [•]
, tr (•) and | • | are respectively the expectation, the trace and the determinant operator. The superscript H denotes the complex conjugate transposition. p is the dimension of the target scattering vector x (p = 3 for the reciprocal case).

The product model

For textured areas, the observed covariance matrix Z can be written as the product of a texture parameter µ with the covariance matrix for homogeneous surface, Z = µZ h . Generally the texture is polarimetric dependent and µ is represented by a matrix. Here, µ is assumed to be a positive scalar parameter. The probability density function of the covariance matrix Z can be derived by the following equation [START_REF] Joughin | Polarimetric Density Functions for Multilook Polarimetric Signatures[END_REF] [3] [START_REF] Lopès | Optimal Speckle Reduction for the Product Model in Multilook Polarimetric SAR Imagery and the Wishart Distribution[END_REF]:

p Z (Z|Σ h , α) = ∞ 0 p Z h (Z|µΣ h ) p µ (µ|α) dµ ( 2
)
where α is the parameter set of texture distribution.

For forest areas, the Gamma distribution has been widely used to model the texture. The associated covariance matrix follows a K distribution, so named because it uses the modified K-Bessel function of the second kind [START_REF] Lee | K-Distribution for Multi-Look Processed Polarimetric SAR Imagery[END_REF].

With high resolution images of man-made objects, the Fisher distribution has been successfully introduced to model the SAR clutter [START_REF] Tison | A New Statistical Model for Markovian Classification of Urban Areas in High-Resolution SAR Images[END_REF]. In next section, the advantage of this distribution family for texture modeling is examined and the covariance matrix pdf is derived analytically.

SCENE TEXTURE DISTRIBUTION

Fisher distribution

The pdf of the Fisher distribution is defined by 3 parameters as:

p µ (µ) = Γ(L + M) Γ(L)Γ(M) L Mm Lµ Mm L-1 1 + Lµ Mm L+M (3)
with L > 0, M > 0. The Fisher distribution is not related to electromagnetic wave theory. It involves the Mellin transform and second kind statistics.

Benefit of Fisher distribution

Fig. 1 shows Gamma, Inverse Gamma and Fisher pdf in the κ 2 /κ 3 plan introduced by [START_REF] Nicolas | Introduction aux Statistiques de deuxième espèce: applications des Logs-moments et des Logs-cumulants à l'analyse des lois d'images radar[END_REF]. Gamma and Inverse Gamma pdf are respectively represented by the blue and red line.

Fisher pdf cover all the space between the blue and red line.

As Fisher distribution is equal to the Mellin convolution of a Gamma distribution by an inverse Gamma distribution, it can be viewed as a generalization of the Gamma distribution.

KummerU distributed covariance matrix

As the Fisher distribution fit well a large range of SAR clutter [START_REF] Tison | A New Statistical Model for Markovian Classification of Urban Areas in High-Resolution SAR Images[END_REF], it is of great interest to derive the covariance matrix pdf for a Fisher distributed texture. This pdf is shown in Eq. 4, it is obtained by substituting Eq. 1 and Eq. 3 into Eq. 2. All the mathematical details can be found in [START_REF] Bombrun | Fisher Distribution for Texture Modeling of Polarimetric SAR Data[END_REF].

p Z (Z|Σ h , L, M, m) = L Lp |Z| L-p π p(p-1) 2 Γ(L) • • • Γ(L -p + 1)|Σ h | L × Γ(L + M) Γ(L)Γ(M) L Mm Lp Γ (Lp + M) U (a, b, z) (4)
The covariance matrix distribution of a Fisher distributed clutter can be expressed as a function of the confluent hypergeometric function of the second kind KummerU (denoted

U (•, •, •)) with z = L tr Σ h -1 Z L Mm , a = Lp + M and b = 1 + Lp -L.
In the following, such pdf will be named KummerU distribution.

HIERARCHICAL SEGMENTATION

In this paper, the hierarchical segmentation algorithm proposed by Beaulieu and Touzi [START_REF] Beaulieu | Segmentation of Textured Polarimetric SAR Scenes by Likelihood Approximation[END_REF] is adapted to the KummerU distributed covariance matrix. The segmentation process can be divided into three steps:

1. Definition of an initial partition.

2. For each 4-connex segments pair, the stepwise criterion, SC, is computed. Then, we find and merge the two segments which minimize the criterion.

3. Stop if the maximum number of merges is reached, otherwise go to step 2.

Stepwise criterion

The criterion used in the hierarchical algorithm is based on the log-likelihood function. The hierarchical algorithm merges the two adjacent segments S i and S j which minimizes the loss of likelihood. The stepwise criterion (SC i,j ) can be expressed as [START_REF] Beaulieu | Segmentation of Textured Polarimetric SAR Scenes by Likelihood Approximation[END_REF]: where MLL(•) is the segment maximum log-likelihood function.

SC i,j = MLL(S i ) + MLL(S j ) -MLL(S i ∪ S j ) (5)
For the KummerU distribution, the maximum log-likelihood function for segment S is derived from Eq. 4. After removing terms that will be cancelled in the stepwise criterion, the log-likelihood function can be rewritten as

MLL(S) = -nL ln |C h | + n ln Γ( L + M) -n ln Γ( L) -n ln Γ( M) + nLp ln L M m + n ln Γ(Lp + M) + Z k ∈S ln U Lp + M; 1 + Lp -L; L tr C h -1 Z k L M m (6)
L, M and m are respectively the estimated of the Fisher parameters L, M and m by the log-cumulants method [5] [9]. C h is the best likelihood estimate of Σ h for the segment S.

Segmentation of a synthetic image

The hierarchical segmentation algorithm proposed by Beaulieu and Touzi [START_REF] Beaulieu | Segmentation of Textured Polarimetric SAR Scenes by Likelihood Approximation[END_REF] has been implemented for the Wishart, the K and the KummerU criterion (Eq. 6) with a synthetic image. Then, the PolSAR dataset is generated by multiplying the texture image by the homogeneous covariance matrix (speckle). The segmentation algorithm is executed with an initial partition where each segment is a bloc of 10 × 10 pixels. The initial partition is composed of 400 segments.

Fig. 3. ROC curves

Tab. 1 shows the Kolmogorov distance between the Fisher distributions used in the simulation. Note that this distance is rather small (from 0.049 to 0.102). Segmentation of the simulated data is a difficult task because the same Wishart distribution is used and texture are quite similar, as can be observed in Fig. 2(b)

Table 1. Kolmogorov distance between the 4 textures. textures 1-2 1-3 1-4 2-3 2-4 3-4 d K
0.049 0.074 0.102 0.063 0.092 0.072

Segmentation result

The evaluation of segmentation result is a difficult problem. Fig. 3 presents curves similar to Receiver Operator Characteristic (ROC) curves. ROC curves plot the detection probability p d versus the false alarm probability p f a . At each iteration, a new partition is produced by merging 2 segments. ROC curves are obtained by calculating p d and p f a for each partition. A description of the method to calculate p d and p f a from a partition P can be found in [START_REF] Bombrun | Fisher Distribution for Texture Modeling of Polarimetric SAR Data[END_REF]. In Fig. 3, the ROC curve for the KummerU distribution is better than the curves for the Wishart distribution and K distribution. When p f a = 0.05, the detection probability is much better for the KummerU criterion (0.85) than for the Wishart and K distribution criterion (0.3) and the corresponding partitions, shown in Fig. 4, contain respectively 5, 29 and 11 segments. This shows the advantage of using a criterion that can discriminate between a large range of texture distributions.

HIERARCHICAL CLASSIFICATION

From an over-segmented partition P of K segments given by the hierarchical segmentation algorithm, a classified image can be obtained by using a hierarchical clustering algorithm. First of all, only the L largest segments (named clusters) of the over-segmented partition are selected. Then, at each iteration, the two more similar clusters are merged. The criterion Here, there is no neighborhood restraint of 4-connexity. Clusters are merged iteratively until to reach the desired number of classes. Finally, the (K -L) small segments are merged to the nearest class for the Wishart distance sense. The proposed unsupervised classification algorithm is applied on the well known 8-look L-band PolSAR image over the Oberpfaffenhofen test site. First, the hierarchical segmentation algorithm is applied to segment the image into 10 000 segments. Then, the 3 567 largest clusters (≥ 20 pixels) are merged iteratively. Fig. 5 shows a colored composition in the Pauli basis of the classification obtained on the Oberpfaffenhofen image for 100 classes. High resolution details are well preserved in the classification map. For example, in urban areas, streets and buildings are isolated.

CONCLUSION

In this paper, the Fisher distribution has been used to model the texture of PolSAR data. The Fisher pdf provide a better texture modeling than the Gamma pdf. With the product model hypothesis, the covariance matrix distribution for a Fisher textured clutter has been mathematically established and expressed as a function of the confluent hypergeometric function of the second kind (KummerU). Next, this distribution has been implemented in a hierarchical segmentation algorithm by calculating the log-likelihood function. Comparisons with the Wishart and K distributions have been done for synthetic textured scenes. ROC curves reveal that the use of an appropriate texture distribution is useful to segment textured PolSAR images.

Finally, an unsupervised hierarchical classification algorithm has been proposed and applied on ESAR L-band Pol-SAR data.

FisherFig. 1 .

 1 Fig. 1. Gamma, Inverse Gamma and Fisher pdf in the κ 2 /κ 3 plan

Fig. 2 .

 2 Fig. 2. (a) image containing the 4 segments (ground truth) and Fisher parameters used in the simulation (b) 4-areas synthetic texture image (200 × 200)

Fig. 2 (

 2 b) shows the simulated texture image (200 × 200 pixels). It contains four 100 × 100 Fisher realizations generated by the inversion of the cumulative repartition function. The speckle is constructed from a Wishart distribution with the same parameters over the whole image (8-look).

Fig. 4 .Fig. 5 .

 45 Fig. 4. Segmentation for a given false alarm probability p f a = 0.05 for (a) KummerU criterion (5 segments), (b) Wishart criterion (29 segments), (c) K criterion (11 segments)


V Roca

Digital contact/presence tracing in FR/Europe: lessons learned after two years

Decentralized risk analysis requires sharing data L

Step 1: server collects (GAEN) or computes (CLÉA) a black list

Step 2: distribute publicly the black list through a CDN (w/o access control)

Step 3: download the black list and compare with local history

Key question: nature of information in the black list?

Decentralized CT threat and privacy

Decentralized risk analysis (e.g., GAEN) requires a smartphone to know pseudonyms of all COVID+ users

The GAEN black list of "diagnosed pseudos" is public (see https://hal.inria.fr/hal-02899412).

No effective and scalable way to add access control Not GDPR friendly… and creates discrimination risks (see: abuse of LUCA by police to find witnesses 2022-01)

On-the-shelf Google's GAEN must be audited too! o why does the "joint statement on CT" (April 19th, 2020) totally omit this trade-off?

• No longer valid in case of higher death rates o essentially a matter of user acceptance + a few tricks (easy upload, etc.)

o however efficiency is complex question… What about evaluation-by-design?

• Centralized scales better and enables technical optimizations

• Acceptance raises different questions… ask sociologists, not only privacy experts

FR

  status • GAEN: AT, IT, LV, DE, PL, DK, IE, HR, ES, NL, FI, EE, PT, MT, BE, CZ, SI, LT, NO, CY (22, of which 19 are interconnected through EFGS) • Centralized: FR, HU (Warning: map not accurate)Contact… or presence tracing?Presence Tracing: have you been in a location where one or more co-located people were later tested COVID+? è CLÉA (FR), LUCA and CWA Event Registration (DE), CrowdNotifier (Swiss)… Contact Tracing: did you met one or more COVID+ persons sufficiently close, sufficiently long? è ROBERT (FR), GAEN (in future?) CT: keep a 14-days history of contacts {contact pseudo, time, distance} PT: keep a 14-days history of locations {QRcode (includes location pseudo), time} Step 1: remember… Step 2: check risk… decentralized: locally (CLÉA for PT and GAEN for CT) centralized: by polling a server that computes risks (ROBERT for CT) Alice Contact… or presence tracing? (3)

  centralized CT designs L • Problem: Bob exposes his full contact list to the server L Alice uploads her pseudonyms COVID+ Bob here are my contacts, is one of them at risk? è side effect: server knows the social graph But there are good centralized CT designs J • Server only knows an unordered list of pseudonyms at risk, nothing else • Bob does not reveal anything to the server (idem GAEN J) è knows only a list of pseudos at risk (e.g., ROBERT) Alice uploads her contact list (but not her own pseudo, useless) COVID+ Bob is my pseudo at risk? About CT threat and privacy modesl if ((you can trust the institutions of your country) && (the use of the app is optional) && (a trusted 3 rd party audits the whole system)) { /* health authority has more control, no sensitive * data shared publicly, GDPR compliant, safe */ use a centralized CT; } else { /* risk of abuse cannot be mitigated, * GDPR non-compliancy is perhaps secondary *Decentralized risk analysis (e.g., CLÉA) requires a smartphone to know locations at risk (pseudo) The CLÉA black list of "location pseudos" is public It's not sensitive health data J and all locations Bob visited remain local (unless COVID+) black list of location pseudos Centralized PT threat and privacy models (2) Centralized risk analysis LUCA collects name/phone number of participants to an event for the HA to contact them by phone if needed Potentially helpful from an epidemiologic viewpoint… …but risky from a privacy viewpoint.

o•

  GAEN prohibited« if we're faced with any kind of trade-off one of the guiding insights that we've used through this process is that social graph is more sensitive and more privacy risky than infection status data.» G. Hogben (Dir. privacy Android) Oct. 2020: https://www.youtube.com/watch?v=0gqZJXOO9Ko[offset 3:35] In March 2020 there were suspicions it could work (simulations)• One year later [1] proved CT saved 4200 -8700 lives in UK in Sept.-Dec. 2020 o to compare to the 32 500 deaths during period o but ∼28% of regular users within population [1] "The epidemiological impact of the NHS COVID-19 app", Nature, Vol 594, 17 June 2021. https://www.nature.com/articles/s41586-021-03606-z Epidemiologic: In FR it's hard to conclude • Missing a sound epidemiologic study! o Usage analytics since June 2021 only (required a decree, published Feb. 2021) o The existing efficiency evaluation document is limited in its methodology L o Raises the question of evaluation-by-design… Epidemiologic: no clue for Presence Tracing • No evaluation as far as I know • In FR, after one month, use of PT compromised by generalization of sanitary pass in bars/restaurants/sport centers • Important to know since CT and PT rely on different virus spreading models… Which one is valid? wearing masks or not? Technical: Centralized/ROBERT more efficient • Scalability in peak contamination phases o one of the 3 goals of Digital CT o ROBERT scales better (computation on server + no need to download a black list) days TAC vs. GAEN (CWA DE app) NB: collected Feb. 17-22 nd , 2022 during Omicron wave, keeping only traffic from TAC (resp. CWA) servers

Forward tracing or backward tracing?

	Assume Alice is tested COVID+…
	Q: where's the
	main risk? • Forward: who might have been contaminated by Alice when
	contagious?
	è inform her contacts
	o QRcode scanning is less demanding than BLE

Step 3: share if COVID+… if Alice agreed, upload to server : GAEN for CT: Alice pseudos (only) CLÉA for PT: scanned QR codes (only) ROBERT for CT: contact list (only) Alice upload… COVID+ Contact… or presence tracing? (4) • different epidemiologic assumptions o direct virus transmission (close to a COVID+ person), or o virus remains in a location in the air for some time… • different technological requirements

• Backward: on what occasion/where could Alice be contaminated?

	Contact Tracing
	Presence Tracing
	(does both J)
	è inform people who were there
	è (potentially) identify who is key index, inform him (was he
	asymptomatic?) and his contacts

•

  Being centralized enables symmetric cryptography… prohibited with GAEN o keep a shared key on the server per registered app the threat model WRT institutions versus Google/Apple o the threat model WRT data that decentralized systems share • GAEN acceptable for COVID… But not with a virus with higher death rate o GAEN not being GDPR compliant is an issue… see: www.coronadetective.eu • Result is anyway a trade-off as different stakeholders have different expectations • Digital CT proved to be efficient (UK / GAEN / 28% of usage)

o enables efficient use of 16 bytes of payload for BLE messages è ROBERT BLE messages have (limited) anti-replay + encrypted country code è ROBERT avoids the "single world" syndrome of the pan-European interconnection of GAEN apps (EFGS service): each app downloads all diagnosed keys of all EU countries • There's no single universal solution… Depends on: o

Centralized CT threat and privacy

Centralized risk analysis (e.g., ROBERT) requires the server to know the contacts of COVID+ users No public black list (neither "diagnosed pseudos" nor "at risk" pseudos) GDPR compliant J Yet there are major assumptions regarding the threat model! list of contacts with COVID+ users « 0x5B08 », 10mn « 0xAE20 », 5mn …

Sovereignty

DÉSIRÉ 3 rd way CT: "identify encounters, not devices"

A paradigm shift: from "public device pseudonyms" (GAEN, ROBERT) to "private encounter pseudonyms":

• that change across time and devices • are private to users who met (eavesdroppers powerless)

Alice Bob compute encounter pseudo: "0x20E17" compute encounter pseudo: "0x20E17"

"DESIRE: Leveraging the best of centralized and decentralized contact tracing systems", ACM DTRAP, 2021.

https://hal.inria.fr/hal-03476799/en

And major benefits:

• better privacy protection • highly flexible: centralized or decentralized risk evaluation, at a country own discretion • full interoperability across different deployments (≠ GAEN or ROBERT)

Thank you…

vincent.roca@inria.fr


Bertrand Fontaine
email: bertrand.fontaine@ua.ac.be

Herbert Peremans

Jan Steckel

3D sparse imaging in biosonar scene analysis

Bats can navigate through foliage thanks to their biosonar. To build an image of their environment they emit ultrasonic pulses and analyse the corresponding echo signals which consist of a complex superposition of a set of highly overlapping reflected pulses coming from the targets at their different positions. We propose to use a sparse representation of echoes in overcomplete dictionaries to perform 3D imaging using biosonar signals. We show that even if the echosignal is a superposition of highly overlapping echoes the algorithm can efficiently reconstruct the image only from a single measurement. Those results can be useful from a robot sensor perspective where the robot has to find a free path. It can also help to understand better how bats can efficiently navigate through dense foliage.

I. INTRODUCTION

Navigation through a realistic outdoor environment, i.e one containing trees, bushes and other natural reflectors, while solely relying on sonar would seem to be an impossible task if it were not routinely solved by bats. A bat flying in a natural environment, e.g. a forest, is confronted with large numbers of (randomly) distributed scatterers of various sizes. Yet, it manages to extract the spatial information it needs in order to manoeuvre with great skill through its environment. To do so, the so-called FM-bat emits short Frequency Modulated (FM) calls (e.g. in Fig. 1(a)) and analyse the corresponding echoes (e.g. in Fig. 1(b)).

If we assume that the environment reflects towards the bat a number of glints (strong echoes) much smaller than the number of samples needed to represent the received sound, this 3D image reconstruction task can be formulated as a sparse problem. To solve it, a redundant dictionary is built using shifted version of a simple reflector echo signal as atoms. The shift in time, which corresponds to a certain distance traveled by the soundwaves, allows a reconstruction along the depth axis. It is similar to what is done for seismic signals analysis [START_REF] Mallat | A wavelet tour of signal processing: the sparse way[END_REF].

The atoms of the dictionary are further filtered as they would be by the bat's Head-Related Transfer Function (HRTF) providing spatial cues [START_REF] Aytekin | The bat head-related transfer function reveals binaural cues for sound localization in azimuth and elevation[END_REF], [START_REF] Wotton | Spatially dependent acoustic cues generated by the external ear of the big brown bat, Eptesicus fuscus[END_REF] about the reflector's position. This spatial dependent filtering provides information about the other two degrees of freedom: azimuth and elevation. This approach is similar to the one proposed in [START_REF] Pearlmutter | Monaural source separation using spectral cues[END_REF] for separating sounds from different sources. The HRTF processing is modelled using a phase array which is optimised with respect to a sparse recovery criterion [START_REF] Fuchs | On sparse representations in arbitrary redundant bases[END_REF]. Using a robotic setup, we show how this 3D imaging technique can localise multiple simple objects very efficiently based on a single measurement. We also apply the technique on more complex reflecting objects and show that, even if the results are still good, other problems arise that are not dependent on the algorithm self but e.g. the choice of the bases or the difference in energy between echoes coming from different directions.

II. MEASUREMENT SETUP

A polaroid emitter [START_REF] Biber | The polaroid ultrasonic ranging system[END_REF] is mounted on a pan/tilt platform. The transmit signal sent to the emitter is a downward sweep from 100 khz down to 40kHz with a duration of 1ms at a sampling frequency of 500kSamples/sec. Due to the non-linearities of the polaroid transducer, the first harmonic is also present in the emitted acoustic signal and the duration is a bit longer. The sent signal, which resembles a typical FM-bat call [START_REF] Surlykke | Echolocation behavior of big brown bats, eptesicus fuscus, in the field and the laboratory[END_REF], can be seen in Fig. 1.

The emitter, which has a piston-like directionality, concentrates the emitted energy mainly in the first lobe pointing forward. Therefore, the region in space where signals of sufficient energy can be reflected from is restricted to a truncated cone as shown in Fig. 2. The maximum aperture considered is 25 • in azimuth and elevation.

The receiver consists of a 16 channel (4x4), 2D array (s. Fig. 3) built using Knowles FG23329-PO7 microphones and connected to a custom made 16 channel A/D system sampling the 16 microphone signals simultaneously at 500kSamples/sec. As explained further, the 16 received signals from the array are combined and this combined signal is bandpass filtered between 115kHz and 35kHz, and downsampled to 250kSamples/sec.

III. DEPTH AXIS

As a first approximation, the echo signal received by a bat from a small object can be considered as a linear superposition of echoes, each one a scaled and delayed version of the emitted call. Hence, if it is assumed that the phase of the emitted pulse is not changed upon reflection, a target image R can be entirely characterised by a sum of scaled and time-shifted Dirac impulses

R(t) = A i δ(t -t i ) (1)
where A i denotes the amplitude of the echo reflected by reflector i and t i = 2 * r i /v s with r i denoting the range along the depth axis of reflector i and v s the speed of sound. This can be extended to the case where the depth axis is longer, i.e. when the glints might come from different objects.

The measured echo S(t) can be written as the convolution of the target impulse response and the emitted call C(t) which in our case is a chirp from 100kHz down to 40kHz. Hence, introducing sampled versions of the signals of interest i.e., the column vector S = [S(0.T s ), S(1.T s ), • • • , S(n.T s )] T denoting the received signal with n the total number of samples and T s the sample period, and analogously for the vectors C and R denoting the emitted call and the target image respectively, yields

S = DR (2)
where D is a n × n matrix whose k-th column contains the sampled version of the emitted call C delayed by k samples.

D =                 C(0.Ts) 0 • • • 0 C(1.Ts) C(0.Ts) • • • 0 C(2.Ts) C(1.Ts) • • • 0 . . . . . . . . . C(m.Ts) C((m -1).Ts) • • • 0 0 C(m.Ts) • • • 0 0 0 • • • 0 . . . . . . . . . 0 0 • • • C(0.Ts)                 (3)
For example (see Fig. 4), if the received signal contains two echoes with amplitudes 1 and 0.3 respectively arriving at times k • T s and p • T s , the 'simplest' way to build S from R using D is to set the k-th and the p-th samples of R equal to 1.0 and 0.3 and zero all the other ones. The matrix D is the dictionary to build S from R. In the ideal case, the shifted signal in the base is a perfect copy of the emitted call. In reality, even a finite sphere filters the reflected echo. The question then arises as what should be the signal taken as a basis. Two choices are tried; the emitted call and an echo coming from a wooden sphere of diameter 6cm. In any case, an object has a finite number of facets reflecting towards the direction of emission. Therefore, the number of echoes in S and thus the number of non-zero entries in R will be small with respect to the number of samples n, R can be considered as sparse.

The problem of finding the sparsest solution R of Eq.2 is relaxed, a non negativity constraint is added and noise is taken into account. The final formulation is:

min R 1 2 S -DR 2 2 + γ R 1 , R ≥ 0 (4)
with γ > 0. To solve this convex quadratic optimisation problem, an interior-point method is used [START_REF] Kim | A method for large-scale l1-regularized least squares[END_REF]. The algorithm, whose code can be found at http://www.stanford.edu/ ∼ boyd/ l1 ls/, is fed with D and its transpose as operators. The computational burden can be alleviated significantly by exploiting the fact that the columns of D are shifted versions of each other To illustrate the approach we create a simulated echo composed of 11 equal amplitude echoes each delayed by respectively 20,80,84,144,152,212,242,284,300,360 and 380 µs. The simulated object is thus 6.1cm long and contains point-reflectors spaced by minimum 0.6mm. We can see in Fig. 5 that for high Signal to Noise Ratio (SNR as defined in [START_REF] Simmons | Echo-delay resolution in sonar images of the big brown bat, eptesicus fuscus[END_REF]) the reconstruction is perfect in both amplitude of the glints and their spacing. Even for lower SNR's the spacing estimation remains very good.

To validate the simulations, small wooden spheres with diameter of 12mm and 6mm were attached in pairs on wires mounted one after the other. They were moved apart along the z-axis so that the closest sphere would not shadow the other. The different bases tried were the reflections of single spheres at the position x=0. The results can be seen in Fig. 6.

IV. REFLECTOR LOCALISATION

From spatial hearing theory it is known that upon arrival at the listener's head the received acoustic signals are filtered by the Head Related Transfer Function (HRTF) [START_REF] Blauert | Spatial Hearing[END_REF]. In particular for bats this HRTF induced filtering is considered to provide important spatial cues that help the bat localise its target [START_REF] Aytekin | The bat head-related transfer function reveals binaural cues for sound localization in azimuth and elevation[END_REF], [START_REF] Wotton | Spatially dependent acoustic cues generated by the external ear of the big brown bat, Eptesicus fuscus[END_REF]. We use the filtering induced by the HRTF to estimate the position of the reflector. Assuming that the HRTF filtering is mostly due to the fact that, before reaching the tympanic membrane, sounds bounce against different parts of the outer ear we model this by delaying slightly the individual channels of the microphone array introduced in Section II and scaling each of them differently. Those 32 parameters (16 delays values and 16 scaling factors) are optimised using a genetic algorithm. Although the model could be optimised to match a real bat HRTF, we choose to use a criterion more specific to sparse recovery, i.e. the Fuchs criterion [START_REF] Fuchs | On sparse representations in arbitrary redundant bases[END_REF].

The spectrograms of four filtered echoes are shown in Fig. 7. Due to the different delays and gains applied to the 16 microphone signals by the delay and sum beamformer, additional spectral cues i.e., notches and peaks, are introduced which depend on the spatial position of the reflector. Note that this approach has similarity with the one proposed in [START_REF] Pearlmutter | Monaural source separation using spectral cues[END_REF]. However, in [START_REF] Pearlmutter | Monaural source separation using spectral cues[END_REF], HRTF filtering is used to separate sound sources and not to localise them.

V. 3D IMAGING

The total dictionary D t is built up by taking the union of the dictionaries D each containing shifted atoms filtered by an HRTF corresponding to position i, one dictionary for each possible target position.

D t = |D 1 |...|D i |D k | (5)
To indicate how such an approach would indeed allow to simultaneously determine the target's biosonar image as well as the target's position, four sets of spheres are arranged as shown in Fig. 8(b). The two first sets are horizontal in the frontal plane and the two other ones are vertical in the back plane. The two planes are separated by 20cm. The space between [-20 • , 20 • ] degrees azimuth and elevation is divided in steps of 2.5 • . The depth axis is divided in 200. The energy in the resulting image from each dictionary is taken as the probability to have a glint at this position. The 3D reconstruction is shown in Fig. 8 This assumption should be relaxed if one wants to estimate the position of a complex objects which might not be consider as a set of point reflectors. To do so, pieces of crumpled paper are fixed to a wire and mounted at different positions and configurations in front of the emitter. Most of the results only reconstruct parts of the objects, namely the parts reflecting most energy, as shown in Fig. 9. Indeed, the directivity of the emitter biases the reconstruction towards the ensonified cone and in addition, all of the facets of the objects do not reflect directly at the emitter either. The technique can however be improved if the emitter is allowed to scan across space and/or move around the objects and partial reconstructions are combined. However, a few favorable cases allowed good reconstruction in 3D based on a single measurement only, as shown in Fig. 10. Those measurements were often with objects elongated along the depth axis.

(b) in spherical coordinates

VI. DISCUSSION

We have shown that sparse analysis could be used to estimate the positions of objects in 3D space. The algorithm performs very good as shown by the close similarity between the signal reconstructed from the sparse coefficients and the real measurement. We believe that the problems encountered are not inherent to the technique proposed but are due to some of the simplifying assumptions made. The scheme works very well when an object is composed of a superposition of base signals but performs less well when analysing complex object echoes. This is probably due to the fact that the facets of a complex object do not reflect exactly as spheres or perfect reflectors. This can be improved by choosing a more appropriate base. Possible solutions to address these problems could consist of averaging the features of various possible objects or learning the dictionary [START_REF] Aharon | K-svd: An algorithm for designing overcomplete dictionaries for sparse representation[END_REF], [START_REF] Lewicki | Learning overcomplete representations[END_REF].

The most important problem with the practical implementation of this scheme is the difference in returned energy between different directions. This can be due to the fact that not all facets of a reflecting object are pointing towards the emitter. However, the difference in emitted energy, due to the directionality of the emitter, limited the extent of the reconstructed space based on a single measurement most in our experiments. As we have seen, the sparse reconstruction will favor the strongest glints and will tend to disregard the weaker ones. Moreover, whereas the difference in energy is logarithmic, the sparse representation returns a linear scale solution. In an engineering context, this can be overcome by using a more omnidirectional emitter which would make sure that equal energy is radiated in all directions. Alternatively, moving around emitter and receiver [START_REF] Walker | One tone, two ears, three dimensions: A robotic investigation of pinnae movements used by rhinolophid and hipposiderid bats[END_REF], as the bat does by moving its pinnae and its head [START_REF] Mogdans | The role of pinna movement for the localization of vertical and horizontal wire obstacles in the greater horseshoe bat, rhinolophus ferrumequinum[END_REF], would allow to combine partial reconstructions by focusing on different regions in space consecutively.

We propose that, apart from its engineering interest, studying echolocation in a sparse context could allow to understand better how the brain processes the information contained in the echoes. Indeed, sparse codes have been widely studied in neuroscience [START_REF] Olshausen | Sparse coding of sensory inputs[END_REF] and have proven to be especially efficient for vision [START_REF]Sparse coding with an overcomplete basis set: A strategy emploeyed by v1?[END_REF]. A scheme similar to the approach proposed here has been successfully appied to neural networks [START_REF] Asari | Sparse representations for the cocktail party problem[END_REF]. Moreover, the very low number of spikes produced in the bat auditory system [START_REF] Sanderson | Neural responses to overlapping fm sounds in the inferior colliculus of echolocating bats[END_REF] could make sparse coding a very efficient alternative to rate coding.

Fig. 1 .

 1 Fig. 1. Emitted call in its time-pressure representation (a) and spectrogram (c). Echo from a complex object in its time-pressure representation (b) and spectrogram (d)

Fig. 2 .

 2 Fig. 2. Due to the directivity of the emitter the covered space is a truncated cone.

Fig. 3 .

 3 Fig. 3. The 2D array consisting of 16 microphones.

Fig. 4 .

 4 Fig. 4. The sparse representation R together with the dictionary D allows to reconstruct the received signal S.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Mean reconstructed target images (50 trials) at 4 SNR conditions: (a):SNR=45dB, (b): SNR=35dB, (c): SNR=25dB and (d): SNR=15dB. The true target image consists of 11 echoes (vertical lines) all having the same amplitude.

Fig. 7 .

 7 Fig. 7. Echoes from a sphere positioned at different angles. The phased array and the associated delay and sum beamformer filters the echoes in a spatially dependent way.

Fig. 8 .

 8 Fig. 8. 3D reconstruction of four objects, each made out of three wooden spheres. The big empty circles are centered on the true positions of the wooden spheres. The filled circles indicate the positions found with the sparse representation. Their radius is proportional to their value in the R vector. (a) Frontal view, (b) 3D view. (c) The measurement, (d) the reconstructed signal using the weights given in the sparse representation and (e) the error signal.

Fig. 9 .

 9 Fig. 9. (a) 3D reconstruction of a complex object. The solid line is the approximate position of the object. The circles indicate the positions found with the sparse representation. Their radius is proportional to their value in the R vector. (b) The measurement, (c) the reconstructed signal using the weights given in the sparse representation and (d) the error.

Fig. 10 .

 10 Fig. 10. (a) 3D reconstruction of a complex object. The solid line is the approximate position of the object. The circles indicate the positions found with the sparse representation. Their radius is proportional to their value in the R vector; (b) side view, (c) frontal view and (d) view from above. (e)The measurement, (f) the reconstructed signal using the weights given in the sparse representation and (g) the error.