<figure><img src="image_1.png" /></figure>

HAL
open science

# Moment-Matching Array Processing Technique for diffuse source estimation

Colin Cros, Laurent Ferro-Famil

► To cite this version:

Colin Cros, Laurent Ferro-Famil. Moment-Matching Array Processing Technique for diffuse source estimation. 2025. hal-05420630

HAL Id: hal-05420630

https://hal.science/hal-05420630v1

Preprint submitted on 17 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Moment-Matching Array Processing Technique for diffuse source estimation

Colin Cros¹, Laurent Ferro-Famil¹,²

¹ ISAE Supaero, University of Toulouse, France

² CESBIO, University of Toulouse, France

firstname.lastname@isae-supaero.fr

December 17, 2025

## Abstract

Direction of Arrival (DOA) estimation is a fundamental problem in signal processing. Diffuse sources, whose power density cannot be represented with a single angular coordinate, are usually characterized based on prior assumptions, which associate the source angular density with a specific set of functions. However, these assumptions can lead to significant estimation biases when they are incorrect. This paper introduces the Moment-Matching Estimation Technique (MoMET), a low-complexity method for estimating the mean DOA, spread, and power of a narrow diffuse source without requiring prior knowledge on the source distribution. The unknown source density is characterized by its mean DOA and its first central moments, which are estimated through covariance matching techniques which fit the empirical covariance of the measurements to that modeled from the moments. The MoMET parameterization is robust to incorrect model assumptions, and numerically efficient. The asymptotic bias and covariance of the new estimator are derived and its performance is demonstrated through simulations.

**Key words:** Array processing; Diffuse source; DOA estimation; Central moments; Misspecified Model

## 1 Introduction

Direction-of-arrival (DOA) estimation is central to many array processing applications, including acoustics, radar, wireless communications, underwater sonar, and synthetic aperture radar (SAR) imaging. Classical high-resolution algorithms, such as MUSIC [1] or ESPRIT [2], characterize a small number of point-like sources using line spectrum estimation [3]. When a source consists of a set of numerous radiating elements densely distributed in space [4, 5] this assumption does not hold, and these techniques lead to erroneous results [6]. Several modifications to the MUSIC and ESPRIT algorithms have been proposed to solve the DOA estimation problem in the case of diffuse sources [7, 8, 9, 10]. However, these subspace methods may face severe limitations in the presence of full-rank source signals. Weighted pseudo-subspace fitting (WPSF) [11] has been proposed to overcome this difficulty, but it comes at the cost of higher computational complexity.

Another limitation of these subspace methods is their need for a source distribution model, which is typically assumed to belong to a specific class of functions. When the shape of a diffuse source density remains unknown, setting an *a priori* candidate might induce model misspecifications, and lead to estimation biases, as this will be demonstrated. Model misSpecifications in DOA estimation have been studied from different perspectives, such as, for instance, the
incorrect modelling of array manifolds [12, 13], or misspecified noise statics [14, 15, 16, 17]. A second-order approximation of the covariance matrix has been proposed to create a low-complexity estimation algorithm [18], though it still uses a model to estimate the spread. The use of higher orders has also been explored: in [19], an iterative algorithm is proposed to identify several diffuse sources based on a Taylor expansion of their distributions and estimation of their non-central moments. However the estimation of the source spread still relies on the use of an a priori shape. Finally, methods that do not rely on assumptions about the source distribution have been proposed for monopulse radar [20, 21, 22], although those methods do use prior information about the source's position.

This paper introduces a Moment-Matching Estimation Technique (MoMET) to estimate three relevant characteristics of a narrow diffuse source: its mean DOA, spread, and power. The MoMET algorithm does not require any prior information about the source distribution. Instead, the diffuse source density is parameterized by its moments.

The main contributions are the following:

1. A model for the covariance matrix of a diffuse source based on the central moments of the source's distribution is introduced. This model does not require knowledge of the source power density's shape, and is linear in all parameters except the mean DOA.

2. A low-complexity estimation algorithm, called MoMET, is derived by fitting the proposed model to the sample covariance matrix. This algorithm can be applied to both uniform linear arrays (ULAs) and non-ULAs.

3. The performance of MoMET is derived in terms of bias and variance, and is validated numerically. Specifically, increasing the number of estimated moments is proved to improve the estimation performance.

This work follows on from a preliminary conference communication by the authors [23]. However, the results presented here are more detailed, self-contained,

and contain new materials; for example; no performance studies or problem analyses were conducted in [23]. While the development of MoMET was motivated by forest SAR tomography applications [24] in which forest density of reflectivity does not follows any simple model, MoMET can be applied to any DOA estimation problem.

The rest of the paper is organized as follows. Section 2 introduces and discusses the estimation problem, then Section 3 presents the new estimation method. Section 4 examines the theoretical performance, and Section 5 assesses this performance with simulations.

**Notation.** In the sequel, lowercase boldface letters denote vectors e.g., $\boldsymbol{x} \in \mathbb{C}^M$, and uppercase boldface variables represent matrices e.g., $\mathbf{A} \in \mathbb{C}^{M \times N}$. The trace, inverse, transpose, Hermitian, and the $(k,l)$ entry of a matrix $\mathbf{A}$ are denoted by tr $\mathbf{A}$, $\mathbf{A}^{-1}$, $\mathbf{A}^{\top}$, $\mathbf{A}^{\mathrm{H}}$, and $[\mathbf{A}]_{k,l}$, respectively. $I_M$ is the $(M \times M)$ identity matrix. The Hadamard (element-wise) product and the Kronecker product of two matrices are denoted by $\mathbf{A} \odot \mathbf{B}$ and $\mathbf{A} \otimes \mathbf{B}$, respectively. For a matrix $\mathbf{A} \in \mathbb{C}^{M \times N}$, $\mathrm{vec} \mathbf{A} \in \mathbb{C}^{MN}$ is the vector obtained by stacking the columns of $\mathbf{A}$. For a vector $\boldsymbol{a} \in \mathbb{C}^M$, diag $\boldsymbol{a} \in \mathbb{C}^{M \times M}$ is the diagonal matrix induced by $\boldsymbol{a}$. The Loewner partial order on positive semi-definite matrices is denoted by $\preceq$. The notation $\mathrm{E}[\cdot]$ denotes the expected value of a random variable and $\|\cdot\|$ the Euclidean norm.

## 2 Problem formulation

Consider a 2D space containing a linear array consisting of $M$ sensors, measuring the signal emitted by a distributed source, covering a narrow spectral band. The signals received by the $M$ elements of the array are gathered into a vector, denoted by $\boldsymbol{x}(t) \in \mathbb{C}^M$, modeled as:

$$ \boldsymbol{x}(t) = \int \boldsymbol{a}(\omega)s(\omega, t)d\omega + \boldsymbol{\epsilon}(t). \quad (1) $$

In (1), $t$ stands for the acquisition index, $s(\omega, t)$ denotes the source signal density at angular position $\omega$,
$\epsilon(t)$ stands for the acquisition noise, and $\mathbf{a}(\omega)$ represents the steering vector, which models the relationship between the DOA and the phase shifts measured by the array elements,

$$ \mathbf{a}(\omega) \triangleq (1 \ e^{j u_2 \omega} \ \dots \ e^{j u_M \omega})^{\mathrm{T}}, \quad (2) $$

where $u_k$ denotes the spacing between the first and the k-th sensor, normalized by half the carrier wavelength.

The expression given in (1) is general and can model signals acquired using any array configuration. In the case of an ULA, with $u_k > u_{k-1}$, the array resolution, $\delta_\omega$, and ambiguous domains, $\Delta\omega_a$, are given by:

$$ \delta_\omega \triangleq \frac{2\pi}{Mdu}, \quad \Delta\omega_a \triangleq \min\left(2\pi, \frac{2\pi}{du}\right), $$

respectively, with $du \triangleq |u_k - u_{k-1}|$. It is assumed in the following that the array is operated in a conventional way, i.e., the inter-element spacing and an appropriate demodulation step, required in case of squinted acquisition geometry, guarantee the measured source density is contained within the array's unambiguous angular domain, i.e., $s(\omega, t) = 0$, $\forall |\omega| > \frac{\Delta\omega_a}{2}$. This hypothesis also applies to non-ULA configurations, using the appropriate resolution and ambiguity definitions provided, e.g., in [25]. For the sake of clarity and without loss of generality, it is assumed in the following that $du \approx 1$ (and $du = 1$ when ULAs are considered).

The noise $\epsilon$ is assumed to follow a zero-mean complex circular normal distribution that is spatially and temporally white, thus $\forall (t, t')$:

$$ E[\epsilon(t)\epsilon(t')^{\mathrm{T}}] = 0, \quad E[\epsilon(t)\epsilon(t')^{H}] = \sigma_{\epsilon}^{2}\delta(t - t')\mathbf{I}_{M}, \quad (3) $$

where $\delta$ denotes the Dirac delta function. The source signal is assumed to be a zero-mean, stationary random process following a normal distribution. Its autocorrelation kernel is characterized by its power density $f$ such that $\forall (\omega, \omega', t, t')$:

$$ E[s(\omega, t)s(\omega', t')] = 0, \quad (4a) $$

$$ E[s(\omega, t)s^*(\omega', t')] = \delta(t - t')\delta(\omega - \omega')f(\omega). \quad (4b) $$

Therefore, the measurements $\mathbf{x}(t)$ follow an unconditional model [26]: they are i.i.d., zero-mean. Their covariance matrix is:

$$ \mathbf{R} \triangleq E[\mathbf{x}(t)\mathbf{x}(t)^{\mathrm{T}}] = \int \mathbf{a}(\omega)\mathbf{a}(\omega)^{\mathrm{H}}f(\omega) \, d\omega + \sigma_{\epsilon}^{2}\mathbf{I}_{M}. \quad (5) $$

The objective of this study is to estimate, without any prior information on the source power density, three characteristics of the source: (i) its power $P$, (ii) its mean value $\omega_0$, and (iii) its spread $\sigma_\omega$. They correspond to the first three moments of $f$:

$$ P \triangleq \int f(\omega) \, d\omega, \quad (6a) $$

$$ \omega_0 \triangleq \frac{1}{P} \int \omega f(\omega) \, d\omega, \quad (6b) $$

$$ \sigma_{\omega}^{2} \triangleq \frac{1}{P} \int (\omega - \omega_{0})^{2} f(\omega) \, d\omega. \quad (6c) $$

The main difficulty is that no model for the power density $f$ is available.

# 3 Moment-Matching Estimation Technique

In order to retrieve the characteristics of the power density without any prior on the shape of the density, this section introduces a new model based on a parameterization of the power density using its central moments. The next subsection provides the motivations behind the MoMET.

## 3.1 Motivation

Estimating a power density, $f$, known to belong to a certain class of functions parameterized by some vector $\psi$, is a classical, well-studied unconditional DOA estimation problem, see e.g., [26], whose ML solution can be derived [27]. However, in the case of an unknown density, a more generic approach is required.

As specified in (5), the entries of $\mathbf{R}$ do not directly observe the power density $f$, but its spectrum. Indeed, the $(k, l)$ entry of $\mathbf{R}$,

$$ [\mathbf{R}]_{k,l} = \int f(\omega)e^{j(u_k - u_l)\omega} \, d\omega + \delta_{k,l}\sigma_{\epsilon}^2, \quad (7) $$
corresponds to the spectrum sampled at the spectral coordinate {$u_k - u_l$}. Therefore, for a limited number of array elements, only the low-pass version of $f$ can be reconstructed from $\mathbf{R}$. Figure 1 presents the low-pass reconstructions obtained with three commonly used power densities.

Using a Dirac function as a very coarse approximation of a narrow source power density leads to a simple covariance model,

$$ \hat{\mathbf{R}}_{\text{Dirac}}(\omega_0, P, \sigma_{\epsilon}^2) = P\mathbf{a}(\omega_0)\mathbf{a}(\omega_0)^H + \sigma_{\epsilon}^2 \mathbf{I}_M. \quad (8) $$

The model (8) can be employed to estimate $\omega_0$, $P$ and $\sigma_{\epsilon}^2$ [28, 29], but is not refined enough to measure the source spread.

In the frequency domain, approximating the power density $f$ by a Dirac impulsion is equivalent to seek for a constant spectrum (up to a spatial shift). As $f$ is narrow, its Fourier transform is indeed “flat” near the origin, nonetheless a constant approximation may be too crude. To allow for finer approximations and capture more complex behaviors, the proposed MoMET approximates the spectrum of $f$ by a polynomial instead of a constant. Furthermore, the Fourier transform of a distribution is a characteristic function, which encodes the distribution’ moments in its Taylor expansion. For small spreads, the characteristic function is very smooth near the origin and some of its descriptors, such as the spread of the source density $\sigma_{\omega}$, can be estimated using a few of its first moments. As an illustration, Figure 2 presents the reconstruction of the power density from the off-diagonal entries of a covariance matrix $\mathbf{R}$ (since the diagonal entries do not contain any information) for different regression orders. A Dirac approximation corresponds to a 0-th order regression. As expected, higher-order regressions produce better reconstructions and the regressed spectra tend to match the true spectrum near the origin, motivating the use of a polynomial model for the spectrum of the power density $f$.

## 3.2 Covariance matrix model

This section introduces a model, $\hat{\mathbf{R}}$, for the covariance matrix, $\mathbf{R}$, based on the central moments of the

power density $f$ which does not assume any prior information on the shape of the power density.

The standard distribution, $p$, is defined as a standardized version of the power density, $f(\omega)$, as:

$$ p(x) \triangleq \frac{\sigma_{\omega}}{P} f(\sigma_{\omega}x + \omega_0) \Leftrightarrow f(\omega) = \frac{P}{\sigma_{\omega}} p\left(\frac{\omega - \omega_0}{\sigma_{\omega}}\right). \quad (9) $$

$p$ is a standard distribution in the sense that it is positive, has unit-mass, zero-mean, and unit-variance:

$$ \int p(x) \, dx = 1, \quad \int xp(x) \, dx = 0, \quad \int x^2 p(x) \, dx = 1. \quad (10) $$

After injecting (9) into (5), the covariance matrix $\mathbf{R}$ becomes:

$$ \mathbf{R} = P\mathbf{a}(\omega_0)\mathbf{a}(\omega_0)^H \odot \mathbf{B} + \sigma_{\epsilon}^2 \mathbf{I}_M, \quad (11a) $$

$$ [\mathbf{B}]_{k,l} \triangleq \int p(x) e^{j(u_k - u_l)\sigma_{\omega}x} dx. \quad (11b) $$

In (11), the position of the source, $\omega_0$, is decoupled from the shape of the distribution, encoded in the shaping matrix $\mathbf{B}$. As $p$ is unknown, a model for $\mathbf{B}$ must be set, or equivalently for the characteristic function of $p$ defined as:

$$ \tilde{p}(\xi) \triangleq \int p(x) e^{j\xi x} dx. \quad (12) $$

Note that by definition of $\mathbf{B}$ (11b), $[\mathbf{B}]_{k,l} = \tilde{p}((u_k - u_l)\sigma_{\omega})$. It is well-known that the $d$-th derivative of the characteristic function, denoted by $\tilde{p}^{(d)}$, evaluated at $\xi = 0$, corresponds to the $d$-th moment of the distribution:

$$ \tilde{p}^{(d)}(0) = j^d \mu_d, \quad \mu_d \triangleq \int x^d p(x) dx. \quad (13) $$

Hence, for any $D$, the $D$-th order Taylor expansion at 0 of $\tilde{p}$ is:

$$ \tilde{p}(\xi) = \sum_{d=0}^{D} \frac{j^d}{d!} \mu_d \xi^d + O(\xi^{D+1}). \quad (14) $$

Since the maximum value at which $\tilde{p}$ is evaluated is $\xi_M = \sigma_{\omega}u_M$, and the exponential power series has
<figure><img src="image_2.png" /><figcaption>Figure 1: Comparison of three power densities and their low-pass versions computed from the entries of $\mathbf{R}$. The parameters used are $\omega_0 = 1$ rad, $\sigma_\omega = 2\pi \times 0.05$ rad, $P = 100$ W, and $\sigma_\epsilon^2 = 10$ W. The array is a ULA with $M = 10$ sensors.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Regressions of the spectrum of a uniform power density $f$ and associated reconstructions. The regression points correspond to the off-diagonal entries of $\mathbf{R}$, $\xi_k = \sigma_\omega u_k$. The same numerical values as in Fig. 1 have been used except $\omega_0 = 0$.</figcaption></figure>
an infinite radius of convergence, the remainder of the expansion is negligible for large D.

It is proposed to parameterize $\tilde{p}$ by retaining the first $D$ coefficients of the Taylor expansion. The choice of the expansion order $D$ is discussed in Section 4. According to (10), the first three moments are $\mu_0 = 1$, $\mu_1 = 0$ and $\mu_2 = 1$. The moment-based model $\hat{p}$ is then defined from a vector of new parameters, $\boldsymbol{\mu} = (\mu_3, \dots, \mu_D) \in \mathbb{R}^{D-2}$, as:

$$ \hat{p}(\xi, \boldsymbol{\mu}) \triangleq \sum_{d=0}^{D} \frac{j^d}{d!} \mu_d \xi^d = 1 - \frac{\xi^2}{2} + \sum_{d=3}^{D} \frac{j^d}{d!} \mu_d \xi^d. \quad (15) $$

Injecting (15) into (11) yields the following models for $\mathbf{R}$ and $\mathbf{B}$:

$$ \hat{\mathbf{R}}(\omega_0, P, \sigma_\omega, \boldsymbol{\mu}, \sigma_\epsilon^2) \triangleq P\mathbf{a}(\omega_0)\mathbf{a}(\omega_0)^H \odot \hat{\mathbf{B}}(\sigma_\omega, \boldsymbol{\mu}) + \sigma_\epsilon^2 \mathbf{I}_M, \quad (16a) $$

$$ \hat{\mathbf{B}}(\sigma_\omega, \boldsymbol{\mu}) \triangleq \mathbf{U}^{\circ 0} - \frac{\sigma_\omega^2}{2} \mathbf{U}^{\circ 2} + \sum_{d=3}^{D} \frac{j^d \sigma_\omega^d \mu_d}{d!} \mathbf{U}^{\circ d}, \quad J(\boldsymbol{\theta}) \triangleq \frac{1}{2} \left\| \mathbf{W}^{1/2} (\bar{\mathbf{R}}_N - \hat{\mathbf{R}}(\boldsymbol{\theta})) \mathbf{W}^{1/2} \right\|^2 = \frac{1}{2} \left\| \bar{\mathbf{R}}_N - \hat{\mathbf{R}}(\boldsymbol{\theta}) \right\|_{\mathbf{W}}^2, \quad (19) $$

where $\mathbf{U}^{\circ d}$ denotes the $d$-th Hadamard power of the array geometry matrix $\mathbf{U}$ whose entries are $[\mathbf{U}]_{k,l} \triangleq u_k - u_l$.

The estimation of the parameters in (16) can be simplified by introducing the vector $\boldsymbol{\nu} = (\nu_2, \dots, \nu_D) \in \mathbb{R}^{D-1}$ whose entries are $\nu_d \triangleq P\mu_d\sigma_\omega^d$. As both $P$ and $\sigma_\omega$ are positive, there is a one-to-one correspondence between $(P, \sigma_\omega, \boldsymbol{\mu})$ and $(P, \boldsymbol{\nu})$, so estimating one or the other is equivalent. However, the parameterization $\boldsymbol{\theta} = (\omega_0, P, \boldsymbol{\nu}, \sigma_\epsilon^2)$ is more convenient, as all the variables, except $\omega_0$, enter linearly in the expression of $\hat{\mathbf{R}}$. Let $\boldsymbol{\alpha} \triangleq (P, \boldsymbol{\nu}, \sigma_\epsilon^2)$ denote the set of linear parameters of $\boldsymbol{\theta}$ and introduce the matrices $\{\mathbf{A}_k\}$ associated with these parameters:

$$ \mathbf{A}_k \triangleq \left\{ \begin{array}{ll} \mathbf{11}^\top & \text{for the parameter } P, \\ \frac{j^d}{d!} \mathbf{U}^{\circ d} & \text{for the parameter } \nu_d, \\ \mathbf{I}_M & \text{for the parameter } \sigma_\epsilon^2. \end{array} \right. \quad (17) $$

Then, grouping all the linear terms in a single matrix

$\hat{\mathbf{B}}(\boldsymbol{\alpha})$ leads to the model for the covariance matrix:

$$ \hat{\mathbf{R}}(\boldsymbol{\theta}) \triangleq \mathbf{a}(\omega_0)\mathbf{a}(\omega_0)^H \odot \hat{\mathbf{B}}(\boldsymbol{\alpha}), \quad (18a) $$

$$ \hat{\mathbf{B}}(\boldsymbol{\alpha}) \triangleq \sum_k \alpha_k \mathbf{A}_k = P \mathbf{U}^{\circ 0} + \sum_{d=2}^{D} \nu_d \frac{j^d}{d!} \mathbf{U}^{\circ d} + \sigma_\epsilon^2 \mathbf{I}. \quad (18b) $$

Note that, since the Taylor expansion is truncated, the moment-based model (18) is misspecified by design. The consequences of this misspecification are discussed in Section 4.

### 3.3 Estimation algorithm

The MoMET algorithm aims at retrieving the parameters $\boldsymbol{\theta}$ in (18) using COMET estimators [30], which are based on a generalized least squares method offering a numerically efficient alternative to the ML estimator. Introduce the cost function,

$$ J(\boldsymbol{\theta}) \triangleq \frac{1}{2} \left\| \mathbf{W}^{1/2} (\bar{\mathbf{R}}_N - \hat{\mathbf{R}}(\boldsymbol{\theta})) \mathbf{W}^{1/2} \right\|^2 = \frac{1}{2} \left\| \bar{\mathbf{R}}_N - \hat{\mathbf{R}}(\boldsymbol{\theta}) \right\|_{\mathbf{W}}^2, \quad (19) $$

where $\bar{\mathbf{R}}_N \triangleq \frac{1}{N} \sum_{l=1}^N \mathbf{x}(t)\mathbf{x}(t)^H$ is the sample covariance matrix and $\mathbf{W}$ is an Hermitian weight matrix. The COMET estimator is defined as the minimizer of $J$:

$$ \hat{\boldsymbol{\theta}} = \arg \min_{\boldsymbol{\theta}} J(\boldsymbol{\theta}). \quad (20) $$

The user-defined weight matrix $\mathbf{W}$ may be considered as a processing parameter providing specific properties to the estimation results. In particular, $\mathbf{W} = \mathbf{I}$ yields an ordinary unweighted least squares estimation, while setting $\mathbf{W} = \bar{\mathbf{R}}_N^{-1}$ leads to an asymptotically efficient<fnref n="1" /> solution [30].

As proposed in [30], (20) may be solved by introducing the $M^2 \times (D+1)$ matrix $\mathbf{L}$ whose $k$-th column $[\mathbf{L}]_k \triangleq \text{vec} \mathbf{A}_k$, so that $\mathbf{L}\boldsymbol{\alpha} = \text{vec} \hat{\mathbf{B}}(\boldsymbol{\alpha})$, and the following notations:

$$ \mathbf{\Phi}(\omega_0) \triangleq \text{diag} \mathbf{a}(\omega_0), \quad (21a) $$

$$ \mathbf{\Psi}(\omega_0) \triangleq \mathbf{\Phi}(\omega_0)^H \otimes \mathbf{\Phi}(\omega_0), \quad (21b) $$

$$ \bar{\mathbf{r}}_N \triangleq \text{vec} \bar{\mathbf{R}}_N. \quad (21c) $$

<footnote n="1">¹if the model is well-specified</footnote>
Using these notations and $\theta = (\omega_0, \alpha)$, the cost function, (19), becomes:

$$J(\theta) = \frac{1}{2} (\bar{r}_N - \Psi(\omega_0) L \alpha)^H (W^\top \otimes W) (\bar{r}_N - \Psi(\omega_0) L \alpha). \qquad (22)$$

As pointed out in [30], for a fixed value of $\omega_0$, (22) is a quadratic form in $\alpha$ whose minimum is achieved at:

$$\hat{\alpha}(\omega_0) \triangleq [(\Psi(\omega_0)L)^H (W^\top \otimes W) \Psi(\omega_0)L]^{-1} (\Psi(\omega_0)L)^H (W^\top \otimes W) \bar{r}_N \qquad (23)$$

Finally, $\omega_0$ can be estimated using a simple uni-dimensional optimization of the following criterion, obtained by injecting (23) into (22):

$$\hat{\omega}_0 = \arg \max \left[ (\Psi(\omega_0)L)^H (W^\top \otimes W) \bar{r}_N \right]^\top \times [(\Psi(\omega_0)L)^H (W^\top \otimes W) \Psi(\omega_0)L]^{-1} \times (\Psi(\omega_0)L)^H (W^\top \otimes W) \bar{r}_N \qquad (24)$$

The linear parameters $\hat{\alpha}$ are then deduced by evaluating (23) at $\hat{\omega}_0$. The parameters of interest are finally obtained as:

$$\hat{P} = \hat{\alpha}_0, \qquad \hat{\sigma}_{\omega}^2 = \hat{\alpha}_1/\hat{\alpha}_0. \qquad (25)$$

The complete estimation scheme is summarized in Algorithm 1. The optimization occurs during Step 2 with the estimation of $\omega_0$. A details some implementation aspects that speed up this optimization.

### Algorithm 1 MoMET algorithm.

1. Form the sample covariance matrix $\bar{R}_N$ from the measurements $\{x(t)\}$.

2. Estimate $\hat{\omega}_0$ by solving the one-dimensional optimization problem (24).

3. Deduce the linear parameters $\hat{\alpha}$ from (23) with $\hat{\omega}_0$.

4. Deduce the source characteristics $P$ and $\sigma_\omega$ from (25).

## 4 Performance analysis

The covariance matrix model $\hat{R}(\theta)$ being designed by truncating the Taylor expansion of the characteristic

function $\tilde{p}$, (15), it is, in general, misspecified. Let $\theta^{\text{true}}$ denote the vector of true parameters,

$$\theta^{\text{true}} = (\omega_0^{\text{true}}, P^{\text{true}}, \nu^{\text{true}}, \sigma_{\epsilon}^{\text{true}}).$$

The model's misspecification has two consequences: (i) there may be no parameter that reconstructs the true covariance $R$, i.e., $\forall \theta$, $\hat{R}(\theta) \neq R$ (in particular, $\hat{R}(\theta^{\text{true}}) \neq R$), and (ii) even if there exists $\theta'$ such that $\hat{R}(\theta') = \hat{R}$, the parameters in $\theta'$ shall not be interpreted as the true ones, i.e., $\theta' \neq \theta^{\text{true}}$.

The choice of the expansion order, $D$, impacts the misspecification, and therefore the estimator performance. Increasing the expansion order $D$ generates finer reconstructions of the characteristic function $\tilde{p}$, and therefore of the covariance matrix $R$, and is expected to produce better estimates. However, the number of parameters to be estimated increases with the expansion's order. Slow estimation convergence may be experienced for higher order values. The choice of $D$ relies on a trade-off between the estimator's bias and its variance, which is discussed in the following two subsections. Subsection 4.1 presents the convergence properties of the estimators, and Subsection 4.2 its asymptotic covariance.

For the sake of clarity throughout this section, the number of samples $N$ is emphasized in the notation of the cost function, which is denoted by $J_N(\theta)$. Furthermore, as the number of samples $N$ goes to infinity, $\bar{R}_N$ is a strongly consistent estimator of $R$. Therefore $J_N$ converge almost surely to the asymptotic cost function,

$$J_{\infty}(\theta) \triangleq \frac{1}{2} \left\| \hat{R}(\theta) - R \right\|_{W}^{2}. \qquad (26)$$

Finally, it should be noted that the order $D$ is constrained by the number of sensors $M$. The parameter vector $\theta = (\omega_0, P, \nu, \sigma_\epsilon^2)$, with $\nu = (\nu_2, \dots, \nu_D) \in \mathbb{R}^{d-1}$, consists of $d+2$ parameters to be identified. Moreover, a general linear array with $M$ sensors can observe up to $M(M-1)+1$ parameters. This number is reduced to $2M-1$ for a ULA. Unambiguous identification requires the number of parameters to be estimated does not exceed the number of independent observations [31, 30]. Therefore, the maximal
possible order is given by:

$$D_{\max} = \begin{cases} M(M-1)-1 & \text{for a general LA,} \\ 2M-3 & \text{for a ULA.} \end{cases} \tag{27}$$

## 4.1 Accumulation point and asymptotic bias

This subsection focuses on the convergence of the MoMET estimator. To derive asymptotic properties, it is first assumed that the asymptotic cost function, $J_{\infty}$, has a global minimum on the ambiguity domain. That assumption is discussed at the end of that subsection. This assumption yields the following result.

**Lemma 1.** If $J_{\infty}$ has a global minimum, then the estimator $\hat{\theta}$ converges almost surely to this accumulation point,

$$\hat{\theta} \xrightarrow{\text{a.s.}} \theta_0 \triangleq \arg\min_{\theta} J_{\infty}(\theta). \tag{28}$$

*Proof.* See B.

<figure><img src="image_4.png" /></figure>

As already mentioned, because of the model's misspecification, the accumulation point, $\theta_0$, generally differs from the true parameter, $\theta^{\text{true}}$, i.e., the estimator is biased. Nonetheless, asymptotic upper bounds for the bias $\|\theta_0 - \theta^{\text{true}}\|$ can be derived as follows.

First, note that by construction of $\hat{p}$, (15), and $\hat{R}(\theta)$, (18):

$$J_{\infty}(\theta^{\text{true}}) = O(\sigma_{\omega}^{2(D+1)}).$$

Furthermore, by definition of $\theta_0$, $0 \le J_{\infty}(\theta_0) \le J_{\infty}(\theta^{\text{true}})$. Thus, $J_{\infty}(\theta^{\text{true}}) - J_{\infty}(\theta_0) = O(\sigma_{\omega}^{2(D+1)})$. Since $\nabla J_{\infty}(\theta_0) = 0$, the Taylor expansion of $J_{\infty}$ around $\theta_0$ yields,

$$J_{\infty}(\theta^{\text{true}}) - J_{\infty}(\theta_0) = (\theta_0 - \theta^{\text{true}})^{\top} \mathbf{H}_{\infty}(\theta_0)(\theta_0 - \theta^{\text{true}}) + o(\|\theta_0 - \theta^{\text{true}}\|^2) \tag{29}$$

where $\mathbf{H}_{\infty}(\theta)$ is the Hessian matrix of the cost function $J_{\infty}$ evaluated at $\theta$. Then, the following lemma allows to deduce the asymptotic behavior of the bias.

**Lemma 2.** There exist $\sigma_0 > 0$ and $C > 0$ such that $\forall \sigma_{\omega} < \sigma_0$, $\mathbf{H}_{\infty}(\theta_0) \succeq CI$, where $\theta_0$ is a function of $\sigma_{\omega}$.

*Proof.* See B.

<figure><img src="image_5.png" /></figure>

According to Lemma 2, there exists a constant $C > 0$ such that for $\sigma_{\omega}$ small enough,

$$(\theta_0 - \theta^{\text{true}})^{\top} \mathbf{H}_{\infty}(\theta_0)(\theta_0 - \theta^{\text{true}}) \ge C \|\theta_0 - \theta^{\text{true}}\|^2. \tag{30}$$

Finally, injecting (30) into (29) leads to $\|\theta_0 - \theta^{\text{true}}\|^2 = O(\sigma_{\omega}^{2(D+1)})$, and thus:

$$\|\theta_0 - \theta^{\text{true}}\| = O(\sigma_{\omega}^{D+1}). \tag{31}$$

Thus, the estimation bias decreases as the source's spread decreases, and when the expansion order $D$ increases. However, the significance of this asymptotic results should be weighed against the slowdown in estimation convergence as adding parameters may require more samples to achieve similar performance, as detailed in Subsection 4.2.

Finally, the existence of a unique accumulation point must be discussed. Examples of the asymptotic cost function (26) are represented in Figure 3 for different orders $D$ and two distributions. First, it can be observed that if the order is set to $D = D_{\max}$, there are several minima for which $R$ is perfectly reconstructed, i.e., such that $J_{\infty}(\omega) = 0$. As this may be observed in Figure 3, in the ULA case, the ambiguous minimum is located at $\omega = \omega_0 + \pi$, while in the non-ULA the period of the ambiguity is not well defined [25]. Such an ambiguity is recurrent with COMET method, for example it is also present in the decoupled estimation algorithm presented in [32, 33]. It can be explained as follows. The parameter $\omega_0$ represents the mean of the power density, $f(\omega)$. Since this function is periodic, shifting it by half a period also produces a zero-mean function. However, this second solution can be discarded because it yields non-physical estimates, such as a negative power. Moreover, it can also be observed in Figure 3a that $J_{\infty}$ may have several minima even when $D < D_{\max}$, e.g., in that figure for $D = 3$. That second observation prevents a general study. However, those minima remain in the neighborhood of $\theta_0$ and all the previous development still applies.
<figure><img src="image_6.png" /><figcaption>Figure 3: Loss function $J_{\infty}(\omega, \hat{\alpha}(\omega))$ for different orders. The same numerical values as in Fig. 1 have been used (in particular $\omega_0 = 1$). In (a) the source has an uniform distribution and the array is uniform with 7 sensors and $u_k = k$ ($D_{\max} = 11$). In (b) the source has a Gaussian distribution and the array is non-uniform with 3 sensors ($D_{\max} = 5$): unevenly distributed sensor positions have been generated by adding a random zero-mean Gaussian shift with variance $\sigma^2 = 0.01$ to the ULA coordinates.</figcaption></figure>

## 4.2 Asymptotic covariance

This section studies the behavior in covariance of the MoMET estimator, $\hat{\theta}$, around its accumulation point, $\theta_0$, for large numbers of samples $N$. The approach used in this work follows the one described in [34] for the ML estimation of misspecified models.

As $\hat{\theta}$ converges in distribution to $\theta_0$, a Taylor expansion of the gradient of $J_N$ around $\theta_0$ yields:

$$ \mathbf{0} = \nabla J_N(\hat{\theta}) = \nabla J_N(\theta_0) + \mathbf{H}_N(\theta_0)(\hat{\theta} - \theta_0) + o\left(\|\hat{\theta} - \theta_0\|\right), \quad (32) $$

where $\mathbf{H}_N(\theta)$ denotes the Hessian of $J_N$ evaluated at $\theta$.

Following the same approach as Cramér [35, pp. 500-503] and White [36, Th. 3.2], in the limit of large samples, (32) gives:

$$ \hat{\theta} - \theta_0 \sim_N \mathbf{H}_\infty(\theta_0)^{-1} \nabla J_N(\theta_0). \quad (33) $$

Then, the asymptotic covariance of $\hat{\theta}$ is:

$$ C_N \triangleq \mathrm{E} \left[ (\hat{\theta} - \theta_0)(\hat{\theta} - \theta_0)^\top \right] \sim_N \frac{1}{N} \mathbf{H}_\infty(\theta_0)^{-1} C_0 \mathbf{H}_\infty(\hat{\theta}, \theta_0) \quad (34) $$

where $C_0 \triangleq \lim_{N \to \infty} NE[\nabla J_N(\theta_0)\nabla J_N(\theta_0)^\top]$.

The details of the computation of $\mathbf{H}_\infty(\theta_0)$ and $C_0$

are given in C. Their entries are:

$$ [\mathbf{H}_\infty(\theta_0)]_{k,l} = \mathrm{tr} \left\{ \mathbf{W} \frac{\partial^2 \hat{\mathbf{R}}(\theta_0)}{\partial \theta_k \partial \theta_l} \mathbf{W} (\hat{\mathbf{R}}(\theta_0) - \mathbf{R}) \right\} \\ + \mathrm{tr} \left( \mathbf{W} \frac{\partial \hat{\mathbf{R}}(\theta_0)}{\partial \theta_k} \mathbf{W} \frac{\partial \hat{\mathbf{R}}(\theta_0)}{\partial \theta_l} \right), \quad (35a) $$

$$ [C_0]_{k,l} = \mathrm{tr} \left( \mathbf{W} \frac{\partial \hat{\mathbf{R}}(\theta_0)}{\partial \theta_k} \mathbf{W} \mathbf{R} \mathbf{W} \frac{\partial \hat{\mathbf{R}}(\theta_0)}{\partial \theta_l} \mathbf{W} \mathbf{R} \right). \quad (35b) $$

The expression of the asymptotic covariance, (34), echoes the expression of the Huber's Sandwich Bound, also known as the Misspecified Cramér-Rao Bound (MCRB), see e.g., [17] for an overview. Under mild conditions [36, 37, 17], the MCRB is the asymptotic covariance of the Misspecified Maximum Likelihood (MML) estimator. Despite the expected resemblance between the two covariances, which were indeed derived through a similar process, it should be noted that they are computed at two different points. The covariance matrix of $\hat{\theta}$ (34) is computed at accumulation $\theta_0$ which minimizes $J_\infty$, while the MCRB is computed at a different accumulation point which minimizes another cost function, namely the Kullback-Leibler distance between the true and the as-
sumed model [38]. This difference implies that, unlike
in the well-specified case [30], the COMET estimator
is not asymptotically equivalent to the (Misspecified)
ML estimator. In particular, they do not converge
toward the same accumulation point since they min-
imize two different metrics.

# 5 Numerical examples

To illustrate the benefits of the MoMET, its perfor-
mance is assessed in simulations and compared with
a parametric ML estimator.

Three distributions, having a Gaussian, uniform,
or exponential shape have been simulated. Unless
stated otherwise, the true parameters are set to
ω₀ = 1 rad, σω = 2π × 0.05 ≈ 0.31 rad, P = 100 W
and σε² = 10 W, and the simulated array is a ULA
with M = 7 (the Fourier resolution of the array
is δω = 2π/M ≈ 0.90 rad). The MoMET algorithm
is run with different orders, D, ranging from 2 to
Dmax − 1 = 10. Two other estimators assuming a
Gaussian power density are also applied for compari-
son: (i) the ML estimator, and (ii) the decoupled es-
timator proposed in [32, 33]. Both are well-specified
and efficient for the Gaussian case, but are misspec-
ified for the uniform and the exponential distribu-
tions. The ML estimator is expected to provide the
best performance but can hardly be applied in real
applications because of its high sensitivity to initial-
ization settings. The decoupled estimator is a numer-
ically stable alternative to the ML; that competitor
has been selected because it has a similar complexity
as the MoMET.

Figure 4 presents the evolution of the asymptotic
bias for different orders as a function of the source's
spread σω. For the MoMET, the biases reduce as
the spread σω reduces, and more precisely they are
O(σω<sup>D+1</sup>), as stated in Section 4.1, (31). As expected,
larger orders result in better asymptotic performance.
For the Gaussian distribution and the uniform distri-
bution, the biases on ω₀ reach the system precision
(~10<sup>-8</sup> for this parameter) which suggests that these
estimators are unbiased. This may be due to the sym-
metry of these distributions.

To assess the practical performance of the estima-

tors, the Root Mean Square Error (RMSE) has been
computed over 1 000 Monte Carlo runs for differ-
ent numbers of snapshots N. Figure 5 compares the
performance of the MoMET applied with different
orders, D. That figure also represents the asymp-
totic covariance and the asymptotic biases. It can be
observed that larger order, D, present larger covari-
ances. However, as they also present smaller bias,
they reach eventually, for large number of samples,
better RMSE performance. This illustrates the well-
known trade-off in estimation between bias and vari-
ance [39]. Figures 6-8 compare the performance of the
MoMET estimator with the competitors on the esti-
mation of the three parameters of interest. Each fig-
ure presents the result for one of the three considered
distributions and three different source's spreads. For
clarity, only the MoMET with minimal order, D = 2,
and maximal order, D = Dmax − 1 = 10, are shown in
the figures. When the source's spread is small com-
pared with the Fourier resolution (e.g., σω = δω/10),
small order MoMET performs as well as the decou-
pled estimation. In this regime, the actual shape of
the distribution does not impact the estimation per-
formance. For larger spreads, applying MoMET with
higher orders produce better results, indicating that
a second-order approximation is no longer sufficient.
The choice of the order should therefore be influenced
by the source's spread. Of course, the competitors
perform better when their distribution assumption is
correct, i.e., for the Gaussian case, Fig. 6. However,
for the two misspecified cases, Figs. 7 and 8, MoMET
reaches similar or better performance than the com-
petitors, especially for the estimation of the spread.

6 Conclusion

This paper has presented a new estimation scheme,
called MoMET, for the characterization of a dif-
fuse source. In order to avoid incorrect models,
the MoMET scheme does not assume any particular
shape for the source’s distribution. Instead, the un-
known distribution is characterized by its mean DOA
and its first central moments. The estimation is then
realized by applying the COMET procedure. A per-
<figure><img src="image_7.png" /><figcaption>Figure 4: Normalized asymptotic biases as a function of the spread $\sigma_{\omega}$ for three shaping distributions. The accumulation points $\theta_0$ were obtained by minimizing $J_{\infty}(\theta)$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5: Normalized RMSE computed over 1 000 Monte Carlo runs for different orders. The dashed lines represent the asymptotic covariances, the dotted lines the asymptotic biases.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Normalized RMSE computed over 1 000 Monte Carlo runs for a source having a Gaussian density. The competitors are well-specified.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7: Normalized RMSE computed over 1 000 Monte Carlo runs for a source having a uniform density. The competitors are misspecified.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8: Normalized RMSE computed over 1 000 Monte Carlo runs for a source having an exponential density. The competitors are misspecified.</figcaption></figure>
formance analysis has derived the asymptotic biases and asymptotic covariance of the MoMET estimator. It has been highlighted that estimating more moments results in a global gain of performance, at the cost of a slower convergence rate. This performance gain has been validated through simulations. The proposed algorithm achieved similar or better performance than a misspecified ML estimator. This work demonstrates how sometimes, using an unrealistic and misspecified but generic model can be better than having an incorrect realistic model.

The proposed model is easy to implement and to optimize: the cost function is a quadratic and all parameters, but one, enter linearly in the expression of the solution. Therefore, the proposed algorithm has a low numerical complexity. Moreover, it can be applied with both uniform or non-uniform linear arrays. All those points make it particularity suitable for large scale imagery applications such as forest SAR tomography.

Several avenues are possible for the future. First, the extension of the MoMET to several sources shall be studied. On a theoretical level, questions remain open and will require further investigation: for example, the absence of bias in the case of symmetric distributions, or the choice of the hyper-parameter. Finally, the upcoming arrival of data from the Biomass mission [40] will enable the proposed algorithms to be applied and validated on real data.

# A Implementation aspects

This appendix discusses numerical aspects of the estimation algorithm, and specifically, the estimation of $\omega_0$ in (20), which represent the most demanding part of the optimization algorithm. In practice, the algorithm is applied a large number of times with different data. The following paragraphs propose to speed up the maximization by reformulating and precalculating some terms.

The optimization problem described in (24) may be reformulated as $\hat{\omega}_0 = \arg \max \mathbf{y}(\omega)^H \mathbf{Y}(\omega)^{-1}$ with:

$$ \mathbf{Y}(\omega) \triangleq (\mathbf{\Psi}(\omega)\mathbf{L})^H (\mathbf{W}^\top \otimes \mathbf{W}) \mathbf{\Psi}(\omega)\mathbf{L}, \quad (36a) $$

$$ \mathbf{y}(\omega) \triangleq \mathbf{L}^H \mathbf{\Psi}(\omega)^H (\mathbf{W}^\top \otimes \mathbf{W}) \bar{\mathbf{r}}_N, \quad (36b) $$

where $\mathbf{\Psi}(\omega)$ is defined in (21), $\mathbf{L}$ is the matrix of linear parameters introduced in Section 3.3, and $\mathbf{W}$ is an Hermitian weighting matrix.

**On the computation of $\mathbf{y}(\omega)$** First, note that for any matrix $\mathbf{Z}$:

$$ [\mathbf{L}^H \mathrm{vec}(\mathbf{Z})]_k = \mathrm{tr} \, \mathbf{A}_k \mathbf{Z}, $$

with $\mathbf{A}_k$ explicitly in (17). In particular as:

$$ \begin{aligned} \mathbf{y}(\omega) &= \mathbf{L}^H \mathbf{\Psi}(\omega)^H (\mathbf{W}^\top \otimes \mathbf{W}) \bar{\mathbf{r}}_N, \\ &= \mathbf{L}^H \mathrm{vec} \{ \mathbf{\Phi}(\omega)^H \mathbf{W} \bar{\mathbf{R}}_N \mathbf{W} \mathbf{\Phi}(\omega) \}, \end{aligned} $$

the k-th entry of $\mathbf{y}(\omega)$ is:

$$ \begin{aligned} [\mathbf{y}(\omega)]_k &= \mathrm{tr} \, \mathbf{A}_k \mathbf{\Phi}(\omega)^H \mathbf{W} \bar{\mathbf{R}}_N \mathbf{W} \mathbf{\Phi}(\omega), \\ &= \mathbf{a}(\omega)^H (\mathbf{W} \bar{\mathbf{R}}_N \mathbf{W} \otimes \mathbf{A}_k^\top) \mathbf{a}(\omega). \end{aligned} $$

The matrices $(\mathbf{W} \bar{\mathbf{R}}_N \mathbf{W} \otimes \mathbf{A}_k^\top)$ are independent from $\omega$, and can be precalculated to speed up the optimization of $\hat{\omega}_0$.

Furthermore, in the special case of a ULA, the computation can be further simplified, as:

$$ \begin{aligned} [\mathbf{y}(\omega)]_k &= \mathbf{a}(\omega)^H (\mathbf{W}_R \otimes \mathbf{A}_k^\top) \mathbf{a}(\omega), \\ &= \sum_{m,n} [\mathbf{A}_k]_{n,m} [\mathbf{W}_R]_{m,n} e^{j(n-m)\omega}, \end{aligned} $$

where $\mathbf{W}_R \triangleq \mathbf{W} \bar{\mathbf{R}}_N \mathbf{W}$. Then, by grouping the terms along the diagonals, and noting that the entries of the $\{\mathbf{A}_k\}$ are constant along the diagonals, one gets:

$$ [\mathbf{y}(\omega)]_k = [\mathbf{A}_k]_{1,1} \mathrm{tr} \, \mathbf{W}_R + 2 \mathrm{Re} \left\{ \sum_{l=1}^{M-1} [\mathbf{A}_k]_{l+1,1} \left( \sum_{m=1}^{M-l} [\mathbf{W}_R]_{m,m+l} \right) e^{jl\omega} \right\}. $$

Thus, there exists a precomputable matrix $\mathbf{B}$ such that:

$$ \mathbf{y}(\omega) = \mathrm{Re} \{ \mathbf{B} \mathbf{a}(\omega) \}. $$

Let $\mathbf{w} = (w_0 \cdots w_{M-1})^\top \in \mathbb{C}^M$ be the vector containing the sum of the diagonals of $\mathbf{W}_R$. The matrix $\mathbf{B}$ is defined as:

$$ [\mathbf{B}]_{k,l} = (1 + \delta(l, 0)) [\mathbf{A}_k]_{l,1} w_l. $$
**On the computation of $Y(\omega)$** Recall that by definition:

$$ [Y(\omega)]_{k,l} = \mathrm{vec}(A_k)^H \Psi(\omega)^H (W^\top \otimes W) \Psi(\omega) \mathrm{vec}(A_l), \quad (37) \\ = \mathrm{tr}\{A_k \Phi(\omega)^H W \Phi(\omega) A_l \Phi(\omega)^H W \Phi(\omega)\}. $$

Although $Y(\omega)$ cannot be simplified in the general case, this is possible in common situations. In the special case of unweighted estimation, i.e., $W = I$, $Y(\omega)$ is independent from $\omega$ and $\bar{R}$:

$$ [Y]_{k,l} = \mathrm{tr} A_k A_l. $$

Furthermore, in the case of a ULA it can be shown that there exists a tensor $C = ((c_{k,l,m})$ such that, for all $k, l$:

$$ [Y(\omega)]_{k,l} = \sum_{m=-2M+2}^{2M-2} c_{k,l,m} e^{jm\omega}. $$

The entries $\{c_{k,l,m}\}$ depend only on $W$ and the $\{A_k\}$ and can also be precomputed.

## B Proof of the lemmas

**Lemma.** If $J_\infty$ has a global minimum, then the estimator $\hat{\theta}$ converges almost surely to this accumulation point,

$$ \hat{\theta} \xrightarrow{\text{a.s.}} \theta_0 \triangleq \arg\min_{\theta} J_\infty(\theta). $$

*Proof.* Let us prove that $\lim_{N \to \infty} \hat{\theta} = \theta_0$ with probability 1. For the sake of clarity, let us emphasize on the number of samples in the definition of $\hat{\theta}$: $\hat{\theta}_N = \arg\min_{\theta} J_N(\theta)$.

Let $\varepsilon > 0$, consider an outcome and let us prove that there exists $N_0$ such that $\forall N \ge N_0$, $\|\theta_0 - \hat{\theta}_N\| \le \varepsilon$. First, note that the search space is a compact set; indeed $\hat{\theta}$ is parameterized only by $\hat{\omega}_0$ which belongs to the ambiguity domain, e.g., $\hat{\omega}_0 \in [-\pi, \pi]$ for a ULA.

Since $J_\infty$ has a unique minimum, $\theta_0$, and $\hat{R}$ is a continuous function of $\theta$, there exists $\eta > 0$ such that, $\forall \theta$:

$$ |J_\infty(\theta) - J_\infty(\theta_0)| \le \eta \implies \|\theta - \theta_0\| \le \varepsilon. $$

Moreover, note that for all $\theta$:

$$ J_N(\theta) = J_\infty(\theta) + \frac{1}{2} \|R - \bar{R}_N\|_W^2 + \langle \bar{R}_N - R, R - \hat{R}(\theta) \rangle_W, \quad (37) $$

where $\langle \cdot, \cdot \rangle_W$ denotes the scalar product induced by $W$. As $\hat{R}(\theta)$ is continuous, $\theta$ belongs to a compact set, $R - \hat{R}(\theta)$ is bounded. Furthermore, as $\bar{R}_N$ converges almost surely toward $R$ according to the Central-Limit Theorem, (37) yields that $J_N$ converges almost surely toward $J_\infty$. Therefore, there exists $N_0$ such that $\forall \theta, \forall N \ge N_0$, $|J_N(\theta) - J_\infty(\theta)| \le \eta/2$. Applying twice that inequality at $\hat{\theta}_N$ and $\theta_0$ gives that for all $N \ge N_0$:

$$ J_\infty(\theta_0) \le J_\infty(\hat{\theta}_N) \le J_N(\hat{\theta}_N) + \frac{\eta}{2} \le J_N(\hat{\theta}_N) + \frac{\eta}{2} \le J_N(\theta_0) + \eta. $$

Therefore, $|J_\infty(\hat{\theta}_N) - J_\infty(\theta_0)| \le \eta$, and thus $\|\hat{\theta}_N - \theta_0\| \le \varepsilon$ as claimed. $\square$

**Lemma.** For almost every linear array, and in particular for a ULA, there exist $\sigma_0 > 0$ and $\varepsilon > 0$ such that $\forall \sigma_\omega < \sigma_0$, $H_\infty(\theta_0) \succeq \varepsilon I$, where $\theta_0$ is implicitly seen as a function of $\sigma_\omega$.

*Proof.* Consider the limit case where $\sigma_\omega = 0$, i.e., the shaping distribution is the Dirac delta function. In that case, the model is not misspecified since $\theta_0 = (\omega_0^{\text{true}}, P^{\text{true}}, 0, \sigma_\epsilon^{2\text{true}}) = \theta^{\text{true}}$, and $R = Pa(\omega_0)a(\omega_0)^H + \sigma_\epsilon^2 I = \hat{R}(\theta_0)$.

In that well-specified case, the entries of the Hessian matrices are (the computation is detailed in C):

$$ [H_\infty(\theta_0)]_{k,l} = \mathrm{tr} \left\{ W \frac{\partial \hat{R}(\theta_0)}{\partial \theta_k} W \frac{\partial \hat{R}(\theta_0)}{\partial \theta_l} \right\}. \quad (38) $$

The derivatives of $\hat{R}(\theta)$ are:

$$ \frac{\partial \hat{R}(\theta_0)}{\partial \theta_k} = a(\omega_0) a(\omega_0)^H \odot A_k, $$

where the matrices $A_k$ are given in (17) for the linear parameters, and $[A_0]_{k,l} = j(u_k - u_l)$ for the parameter $\omega_0$. Therefore, (38) can be reexpressed as:

$$ [H_\infty(\theta_0)]_{k,l} = \mathrm{tr} \{ \Phi(\omega_0)^H W \Phi(\omega_0) A_k \Phi(\omega_0)^H W \Phi(\omega_0) A_l \}, $$
which is the inner product induced by $\Phi(\omega_0)^H W \Phi(\omega_0)$. So, the Hermitian matrix is the Gram matrix associated with the family of matrices $A_k$. This family is linearly independent (it is related with the Vandermonde matrix), so $H_\infty(\theta_0)$ is positive definite. The result is finally obtained by continuity as $\sigma_\omega$ tends to 0.

<figure><img src="image_12.png" /></figure>

### C Computation of the asymptotic covariance

Recall the expression of the asymptotic cost function, (26):

$$J_{\infty}(\theta) = \frac{1}{2} \text{tr} \left\{ W(\hat{R}(\theta) - R)W(\hat{R}(\theta) - R) \right\}.$$

Then, the gradient and the Hessian of $J_\infty$ are:

$$\frac{\partial J_{\infty}(\theta)}{\partial \theta_k} = \text{tr} \left\{ W \frac{\partial \hat{R}(\theta)}{\partial \theta_k} W(\hat{R}(\theta) - R) \right\}, \\ [H_\infty(\theta)]_{k,l} = \text{tr} \left\{ W \frac{\partial^2 \hat{R}(\theta)}{\partial \theta_k \partial \theta_l} W(\hat{R}(\theta) - R) \right\} \\ + \text{tr} \left\{ W \frac{\partial \hat{R}(\theta)}{\partial \theta_k} W \frac{\partial \hat{R}(\theta)}{\partial \theta_l} \right\}.$$

Similarly, the gradient of $J_N$ is:

$$\frac{\partial J_N(\theta)}{\partial \theta_k} = \text{tr} \left\{ W \frac{\partial \hat{R}(\theta)}{\partial \theta_k} W(\hat{R}(\theta) - \bar{R}_N) \right\}, \\ = \text{tr} \left\{ W \frac{\partial \hat{R}(\theta)}{\partial \theta_k} W \hat{R}(\theta) \right\} \\ - \frac{1}{N} \sum_t x(t)^H W \frac{\partial \hat{R}(\theta)}{\partial \theta_k} W x(t)$$

Therefore,

$$\frac{\partial J_N(\theta)}{\partial \theta_k} \frac{\partial J_N(\theta)}{\partial \theta_l} = \text{tr} \left\{ W \frac{\partial \hat{R}(\theta)}{\partial \theta_k} W(\hat{R}(\theta) - \bar{R}_N) \right\} \\ \times \text{tr} \left\{ W \frac{\partial \hat{R}(\theta)}{\partial \theta_l} W \hat{R}(\theta) \right\} \\ - \text{tr} \left\{ W \frac{\partial \hat{R}(\theta)}{\partial \theta_k} W \hat{R}(\theta) \right\} \text{tr} \left\{ W \frac{\partial \hat{R}(\theta)}{\partial \theta_l} W \bar{R}_N \right\} \\ + \frac{1}{N^2} \sum_{t,t'} x(t)^H W \frac{\partial \hat{R}(\theta)}{\partial \theta_k} W x(t) x(t')^H W \frac{\partial \hat{R}(\theta)}{\partial \theta_l} W x(t').$$

Using the fact that the $x(t)$ are i.i.d. and the following property,

$$E\{x(t)^H A x(t) x(t)^H B x(t)\} = \text{tr}\{RA\} \text{tr}\{RB\} + \text{tr}\{RARB\},$$

yields

$$E \left\{ \frac{\partial J_N(\theta)}{\partial \theta_k} \frac{\partial J_N(\theta)}{\partial \theta_l} \right\} = \frac{\partial J_\infty(\theta)}{\partial \theta_k} \frac{\partial J_\infty(\theta)}{\partial \theta_l} \\ + \frac{1}{N} \text{tr} \left( W \frac{\partial \hat{R}(\theta)}{\partial \theta_k} W R W \frac{\partial \hat{R}(\theta)}{\partial \theta_l} W R \right).$$

Finally, since by definition of $\theta_0$, $\nabla J_\infty(\theta_0) = 0$, the first term cancels out and (34) is obtained.

### References

[1] R. Schmidt, Multiple emitter location and signal parameter estimation, IEEE Transactions on Antennas and Propagation 34 (3) (1986) 276–280. doi:10.1109/TAP.1986.1143830.

[2] R. Roy, T. Kailath, Esprit-estimation of signal parameters via rotational invariance techniques, IEEE Transactions on Acoustics, Speech, and Signal Processing 37 (7) (1989) 984–995. doi:10.1109/29.32276.

[3] P. Stoica, R. L. Moses, et al., Spectral analysis of signals, Vol. 452, Pearson Prentice Hall Upper Saddle River, NJ, 2005.
[4] E. Aguilera, M. Nannini, A. Reigber, Wavelet-based compressed sensing for sar tomography of forested areas, IEEE Transactions on Geoscience and Remote Sensing 51 (12) (2013) 5283-5295. doi:10.1109/TGRS.2012.2231081.

[5] Y. Huang, J. Levy-Vehel, L. Ferro-Famil, A. Reigber, Three-dimensional imaging of objects concealed below a forest canopy using sar tomography at l-band and wavelet-based sparse estimation, IEEE Geoscience and Remote Sensing Letters 14 (9) (2017) 1454-1458. doi:10.1109/LGRS.2017.2709839.

[6] T.-P. Jantti, The influence of extended sources on the theoretical performance of the music and esprit methods: narrow-band sources, in: [Proceedings] ICASSP-92: 1992 IEEE International Conference on Acoustics, Speech, and Signal Processing, Vol. 2, 1992, pp. 429-432 vol.2. doi:10.1109/ICASSP.1992.226028.

[7] S. Valaee, B. Champagne, P. Kabal, Parametric localization of distributed sources, IEEE Transactions on Signal Processing 43 (9) (1995) 2144-2153. doi:10.1109/78.414777.

[8] Y. Meng, P. Stoica, K. Wong, Estimation of the directions of arrival of spatially dispersed signals in array processing, IEE Proc-Radar, Sonar Navig 143 (1) (1996).

[9] Q. Wu, K. Wong, Y. Meng, W. Read, Doa estimation of point and scattered sources - veterans, in: IEEE Seventh SP Workshop on Statistical Signal and Array Processing, 1994, pp. 365-368. doi:10.1109/SSAP.1994.572519.

[10] S. Shahbazpanahi, S. Valaee, M. Bastani, Distributed source localization using esprit algorithm, IEEE Transactions on Signal Processing 49 (10) (2001) 2169-2178. doi:10.1109/78.950773.

[11] M. Bengtsson, B. Ottersten, A generalization of weighted subspace fitting to full-rank models, IEEE Transactions on Signal Processing 49 (5) (2001) 1002-1012. doi:10.1109/78.917804.

[12] A. Swindlehurst, T. Kailath, A performance analysis of subspace-based methods in the presence of model errors. i. the music algorithm, IEEE Transactions on Signal Processing 40 (7) (1992) 1758-1774. doi:10.1109/78.143447.

[13] A. Ferréol, P. Larzabal, M. Viberg, Statistical analysis of the music algorithm in the presence of modeling errors, taking into account the resolution probability, IEEE Transactions on Signal Processing 58 (8) (2010) 4156-4166. doi:10.1109/TSP.2010.2049263.

[14] P. Tsakalides, C. Nikias, The robust covariation-based music (roc-music) algorithm for bearing estimation in impulsive noise environments, IEEE Transactions on Signal Processing 44 (7) (1996) 1623-1633. doi:10.1109/78.510611.

[15] T.-H. Liu, J. Mendel, A subspace-based direction finding algorithm using fractional lower order statistics, IEEE Transactions on Signal Processing 49 (8) (2001) 1605-1613. doi:10.1109/78.934131.

[16] Q. Tian, T. Qiu, J. Li, R. Li, Robust adaptive doa estimation method in an impulsive noise environment considering coherently distributed sources, Signal Processing 165 (2019) 343-356. doi:https://doi.org/10.1016/j.sigpro.2019.07.014.

[17] S. Fortunati, F. Gini, M. S. Greco, The unspecified cramer-rao bound and its application to scatter matrix estimation in complex elliptically symmetric distributions, IEEE Transactions on Signal Processing 64 (9) (2016) 2387-2399. doi:10.1109/TSP.2016.2526961.

[18] M. Bengtsson, B. Ottersten, Low-complexity estimators for distributed sources, IEEE Transactions on Signal Processing 48 (8) (2000) 2185-2194. doi:10.1109/78.851999.

[19] S. Shahbazpanahi, S. Valaee, A. Gershman, A covariance fitting approach to parametric localization of multiple incoherently distributed sources, IEEE Transactions on Signal Processing 52 (3) (2004) 592-600. doi:10.1109/TSP.2003.822352.
[20] L. Milstein, Maximum-likelihood estimation of the angular position and extent of a target, IEEE Transactions on Information Theory 27 (2) (1981) 187-199.
doi:10.1109/TIT.1981.1056330.

[21] A. Monakov, Maximum-likelihood estimation of parameters of an extended target in tracking monopulse radars, IEEE Transactions on Aerospace and Electronic Systems 48 (3) (2012) 2653-2665. doi:10.1109/TAES.2012.6237615.

[22] U. Nickel, E. Chaumette, P. Larzabal, Estimation of extended targets using the generalized monopulse estimator: Extension to a mixed target model, IEEE Transactions on Aerospace and Electronic Systems 49 (3) (2013) 2084-2096.
doi:10.1109/TAES.2013.6558043.

[23] C. Cros, L. Ferro-Famil, Moment-based characterization of spatially distributed sources in SAR tomography, arXiv preprint arXiv:2503.12975 (2025).

[24] H. Aghababaei, G. Ferraioli, L. Ferro-Famil, Y. Huang, M. Mariotti D'Alessandro, V. Pascazio, G. Schirinzi, S. Tebaldini, Forest SAR tomography: Principles and applications, IEEE Geoscience and Remote Sensing Magazine 8 (2) (2020) 30-45.
doi:10.1109/MGRS.2019.2963093.

[25] P. Babu, P. Stoica, Spectral analysis of nonuniformly sampled data – a review, Digital Signal Processing 20 (2) (2010) 359-378.
doi:https://doi.org/10.1016/j.dsp.2009.06.019

[26] P. Stoica, A. Nehorai, Performance study of conditional and unconditional direction-of-arrival estimation, IEEE Transactions on Acoustics, Speech, and Signal Processing 38 (10) (1990) 1783-1795.
doi:10.1109/29.60109.

[27] T. Trump, B. Ottersten, Estimation of nominal direction of arrival and angular spread using an array of sensors, Signal Processing 50 (1) (1996) 57-69, special Issue on Subspace Methods, Part I: Array Signal

Processing and Subspace Computations.
doi:https://doi.org/10.1016/0165-1684(96)00003-5.

[28] P. Stoica, K. Sharman, Maximum likelihood methods for direction-of-arrival estimation, IEEE Transactions on Acoustics, Speech, and Signal Processing 38 (7) (1990) 1132-1143.
doi:10.1109/29.57542.

[29] B. Ottersten, M. Viberg, P. Stoica, A. Nehorai, Exact and Large Sample Maximum Likelihood Techniques for Parameter Estimation and Detection in Array Processing, Springer Berlin Heidelberg, Berlin, Heidelberg, 1993, pp. 99-151.
doi:10.1007/978-3-642-77347-1_4.

[30] B. Ottersten, P. Stoica, R. Roy, Covariance matching estimation techniques for array signal processing applications, Digital Signal Processing 8 (3) (1998) 185-210.
doi:https://doi.org/10.1006/dspr.1998.0316.

[31] M. Wax, J. Sheinvald, A. Weiss, Detection and localization in colored noise via generalized least squares, IEEE Transactions on Signal Processing 44 (7) (1996) 1734-1743.
doi:10.1109/78.510620.

[32] O. Besson, P. Stoica, Decoupled estimation of doa and angular spread for a spatially distributed source, IEEE Transactions on Signal Processing 48 (7) (2000) 1872-1882.
doi:10.1109/78.847774.

[33] A. Zoubir, Y. Wang, P. Chargé, A modified comet-exp method for estimating a scattered source, Signal Processing 86 (4) (2006) 733-743.
doi:https://doi.org/10.1016/j.sigpro.2005.07.013.

[34] M. L. Diong, E. Chaumette, F. Vincent, On the efficiency of maximum-likelihood estimators of misspecified models, in: 2017 25th European Signal Processing Conference (EUSIPCO), 2017, pp. 1455-1459.
doi:10.23919/EUSIPCO.2017.8081450.

[35] H. Cramér, Mathematical Methods of Statistics (PMS-9), Princeton University Press, 1999.
[36] H. White, Maximum likelihood estimation of misspecified models, Econometrica 50 (1) (1982) 1–25.

[37] Q. H. Vuong, Cramér-rao bounds for misspecified models, in: Social Sci. Work. Paper, Vol. 652, California Institute of Technology, 1986, pp. 1–38.

[38] H. Akaike, Information theory and an extension of the maximum likelihood principle, in: 2nd International Symposium on Information Theory, 1973, 1973.

[39] A. Hero, J. Fessler, M. Usman, Exploring estimator bias-variance tradeoffs using the uniform cr bound, IEEE Transactions on Signal Processing 44 (8) (1996) 2026–2041. doi:10.1109/78.533723.

[40] S. Quegan, T. Le Toan, J. Chave, J. Dall, J.-F. Exbrayat, D. H. T. Minh, M. Lomas, M. M. D'Alessandro, P. Paillou, K. Papathanassiou, F. Rocca, S. Saatchi, K. Scipal, H. Shugart, T. L. Smallman, M. J. Soja, S. Tebaldini, L. Ulander, L. Villard, M. Williams, The european space agency biomass mission: Measuring forest above-ground biomass from space, Remote Sensing of Environment 227 (2019) 44–60. doi:https://doi.org/10.1016/j.rse.2019.03.032.