Robust Detection for Mills Cross Sonar

Olivier Lerda, Ammar Mian, Guillaume Ginolhac, Jean-Philippe Ovarlez, Didier Charlot

► To cite this version:

Olivier Lerda, Ammar Mian, Guillaume Ginolhac, Jean-Philippe Ovarlez, Didier Charlot. Robust Detection for Mills Cross Sonar. IEEE Journal of Oceanic Engineering, 2024, 49 (3), 10.1109/JOE.2024.3374958. [hal-04651458](https://hal.science/hal-04651458v1)

HAL Id: hal-04651458

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

Submitted on 17 Jul 2024

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Robust Detection for Mills Cross Sonar

Olivier Lerda<fnref n="1" />, Ammar Mian<fnref n="2" />, Guillaume Ginolhac<fnref n="2" />, Senior Member, IEEE, Jean-Philippe Ovarlez<fnref n="3" />,<fnref n="4" />, Member, IEEE, Didier Charlot<fnref n="1" />

**Abstract**—Multi-array systems are widely used in sonar and radar applications. They can improve communication speeds, target discrimination, and imaging. In the case of a multibeam sonar system that can operate two receiving arrays, we derive new adaptive to improve detection capabilities compared to traditional sonar detection approaches. To do so, we more specifically consider correlated arrays, whose covariance matrices are estimated up to scale factors, and an impulsive clutter. In a partially homogeneous environment, the 2-step Generalized Likelihood ratio Test (GLRT) and Rao approach lead to a generalization of the Adaptive Normalized Matched Filter (ANMF) test and an equivalent numerically simpler detector with a well-established texture Constant False Alarm Rate (CFAR) behavior. Performances are discussed and illustrated with theoretical examples, numerous simulations, and insights into experimental data. Results show that these detectors outperform their competitors and have stronger robustness to environmental unknowns.

**Index Terms**—Sonar target detection, Adaptive Normalized Matched Filter, Multiple-Input Multiple-Output, Complex Elliptically Symmetric distributions, Tyler's M-estimator, robustness.

## I. INTRODUCTION

### A. Background and motivations

Forward-Looking sonars are solutions for perceiving the underwater environment. In the context of a growing need for decision-making autonomy and navigation safety, they have become fundamental tools for understanding, anticipating obstacles and potential dangers, analyzing and identifying threats. They offer efficient results, allowing detection, tracking, and classification of surface [1], water column [2], or bottom targets [3], in civil [4] or military applications such as mine warfare [5].

At the detection level, monovariate statistical tests under the Gaussian or non Gaussian interference assumption, defined a priori, remain the prevalent approaches [6]–[8].

Nevertheless, many works on multivariate statistics have shown their great interest compared to algorithms developed from monovariate statistics in a large number of application fields. Indeed, multivariate statistics allow advanced modeling of propagation environments. By following these precepts, [9] gets a central detector with the total unknown of the noise parameters, [10] first derives a detector under the assumption of a known covariance

then substitutes it, through a two-step procedure, by an appropriate estimator, finally, [11] and [12] have shown the relevance of subspace data models in consideration of mismatched signals for which, as an example, the target echo would not come precisely from the center of the main beam. These seminal works are now references in remote sensing, and ground or air surveillance but mainly for radar systems.

Moreover, in the radar field, phenomenal progress has also been made in recent decades, guided by increasingly complex systems which required profound changes in concepts and processing. This is especially the case of Space-Time Adaptive Processing (STAP) for airborne radars [13], which bring considerable improvements in the ability to discriminate moving targets at very low speeds, or Multiple-Input Multiple-Output (MIMO) systems that advance detection performance [14], [15] and resolution [16], [17] by exploiting spatial, frequency or waveform diversities.

Although some preliminary work has emerged in recent years [18]–[24], these methods still seem to be underused in sonar systems.

This paper focuses on the adaptive detection of a point target by a correlated orthogonal arrays sonar system. Inspired by these previous works, we will first show that multibeam systems are perfectly adapted to multivariate formalism. We will then propose two new detectors following the GLRT, and Rao two-step approaches [25]–[27], assuming heterogeneous or partially homogeneous clutter [28], [29].

The performance in a Gaussian environment will first be evaluated. We will show that considering a sonar system with Mills [30] cross arrays leads to a better detectability of targets by reducing the clutter ridge.

Nevertheless, complex multi-normality can sometimes be a poor approximation of physics. This is the case for fitting high-resolution clutter, impulsive noise, outliers, and interference. The Complex Elliptic Symmetric (CES) distributions [31] including the well-known compound Gaussian subclass are then natural extensions allowing the modeling of distributions with heavy or light tails in radar [32]–[35] as in sonar [36]. Mixtures of Scaled Gaussian (MSG) distributions [37] are derived and easily tractable approaches.

In this context, particular covariance matrix estimators are recommended for adaptive processing such as the Tyler estimator [38] or Huber's M-estimator [39]. Their uses lead to very substantial performance gains in [40], [41], and [35]. In our application, these considerations will allow us to design a new covariance matrix estimator. The

<footnote n="1">¹Exail Sonar Systems division, La Ciotat, France.</footnote>

<footnote n="2">²LISTIC, University Savoie Mont-Blanc, Annecy, France.</footnote>

<footnote n="3">³ONERA - DEMR, University Paris-Saclay, Palaiseau, France.</footnote>

<footnote n="4">⁴CentraleSupelec - SONDRA, University Paris-Saclay, Gif-sur-Yvette, France.</footnote>
performance of the detectors in a non-Gaussian impulsive
environment can then be studied. On this occasion, we will
show on experimental data, this estimator's interest in the
robustness to corruption of training data.

B. Paper organization

This paper is organized as follows: Section II presents a
dual array sonar system and the experimental acquisition
conditions on which this work is based. In Section III,
the signal model and detection problem are formalized.
According to the two-step GLRT and Rao test design,
coherent adaptive detectors are derived in Section IV.
The performances are evaluated, compared, and analyzed
in Sections V and VI. Conclusions are given in Section
VII. Proofs and complementary results are provided in
the Appendices.

*Notations:* Matrices are in bold and capital, vectors in bold. Re(.) and Im(.) stand respectively for real and imaginary part operators. For any matrix **A** or vector, **A**^T is the transpose of **A** and **A**^H is the Hermitian transpose of **A**. **I** is the identity matrix and CN(μ, Γ) is the circular complex Normal distribution of mean μ and covariance matrix Γ. ⊗ denotes the Kronecker product.

II. SEAPIX SYSTEM

A. Generalities

The SEAPIX system is a three-dimensional multibeam
echosounder developed by the sonar systems division of
Exail (formerly iXblue) [42]. It is traditionally used by
fishing professionals as a tool to assist in the selection
of catches and the respect of quotas [43], by hydro-
acousticians for the monitoring of stocks and morpho-
logical studies of fish shoals [44], by hydrographers for
the establishment of bathymetric and sedimentary marine
charts [45], [46].

Two uniform linear arrays of 64 elements, arranged in
Mills cross, are entirely symmetric, reversible in transmis-
sion/reception, and electronically steerable. They generate
transverse (i.e. across-track) or longitudinal (along-track)
acoustic swaths of 120° by 1.8°, tiltable on +/-60°, pro-
viding a volumic coverage of the water column.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Multiswath capabilities seen from a schematic representation (left): transverse swath footprint is blue, a 60° steered transverse swath is displayed in red, the longitudinal swath is green, and a 45° steered longitudinal swath is orange. Illustration from the operator software (right): An across-track, an along-track, and a tilted along-track swath are observed, as well as an aggregation of fishes and an already constructed bathymetric map.</figcaption></figure>

B. FLS experiment

The SEAPIX system is being experimented with a
Forward-Looking Sonar configuration for predictive target
detection and identification. In this context of use, the
active face is oriented in the “forward” direction rather
than toward the seabed.

In our study, the sensor is installed on the DriX Un-
crewed Surface Vehicle (USV), and inclined by 20° accord-
ing to the pitch angle to the sea surface (Figure 2 left). In
transmission, the vertical antenna (formerly longitudinal)
generates an enlarged beam of 9° in elevation by 120°
in azimuth. A motion stabilized, and electronically tilted
firing angle allows the upper limit of the -3 dB transmit
beamwidth to graze the sea surface and the lower limit
to reach a 50 m depth bottom at about 300 m range.
In receive, the horizontal antenna (formerly transverse)
generates beams of 2° in azimuth by 120° in elevation and
the vertical antenna (which is used again) of 2° in elevation
by 120° in azimuth. A rigid sphere of 71 cm diameter
(Figure 2 right) is also immersed at 25 m depth in the
middle of the water column.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Exail's DriX USV (left): The cross-shaped linear arrays are visible in the gondola. The real target in open water (right): a metallic sphere of target strength TS = -15 dB.</figcaption></figure>

So after each transmission of 20 ms Linear Frequency
Modulation pulses centered at 150 KHz with a 10 KHz
bandwidth and a Pulse Repetition Interval of 0.5 s, the
sensor signals from the two antennas are simultaneously
recorded, allowing an azimuth and elevation representa-
tion of the acoustic environment (Figure 3).

C. Pre-processing and data format

The signals from the 128 sensors provided by the sonar's
embedded software are demodulated in baseband (In-
Phase and Quadrature components) and decimated at
43 KHz. During reception pre-processing, the digitized
samples are compensated for the time-varying gain, fil-
tered to the bandwidth of the waveform, pulse compressed,
then decimated again to bandwidth.

Finally, a ping dataset is a matrix of about 6000 range
bins from 15 m to 400 m, by 64 sensors, by two arrays.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Azimuth (top) and elevation (bottom) view from a single ping, in 50 m depth in La Ciotat bay.</figcaption></figure>

### III. DETECTION SCHEMES

#### A. Data model for a single array

At each time, we, have acquired two digitalized synchronized complex data vectors of $m = 64$ elements, called “snapshots”, which can be written as (by omitting the temporal parameterization):

$$ \mathbf{x}_i = [x_{i,1} \ x_{i,2} \ \dots \ x_{i,m}]^T \quad (1) $$

where $i = 1, 2$ is the array identifier (respectively horizontal and vertical).

After the pre-processing steps, a point-like target observed on the antenna $i$ is simply:

$$ \mathbf{x}_i = \alpha_i \mathbf{p}_i + \mathbf{z}_i \quad (2) $$

where $\mathbf{x}_i \in \mathbb{C}^m$ is the received signal, $\alpha_i$ is an unknown complex target amplitude, $\mathbf{p}_i \in \mathbb{C}^m$ stands for the known deterministic angular steering vector, $\mathbf{z}_i \in \mathbb{C}^m$ is a mixture of scaled Gaussian (MSG) random vector admitting the stochastic representation:

$$ \mathbf{z}_i \stackrel{d}{=} \sqrt{\tau_i} \mathbf{c}_i. \quad (3) $$

The texture $\tau_i$ is an unknown positive deterministic scalar parameter, presumably different for each range bin (i.e. for all time samples). The speckle $\mathbf{c}_i \sim \mathcal{CN}(\mathbf{0}, \sigma_i^2 \mathbf{M}_{ii}) \in \mathbb{C}^m$ is a complex circular Gaussian random vector whose covariance matrix $\mathbf{M}_{ii}$ is known up to a scaling factor $\sigma_i^2$. The term speckle should be understood in the sense of CES statistics rather than as the result of a sum of contributions from reflections in a resolution cell. This model is strongly related to the class of compound Gaussian distributions [31], which assumes a speckle-independent random texture with a given density $p_\tau$. The MSG distribution is more robust than the Gaussian one because the relative scaling between samples allows flexibility in the presence of heterogeneities, such as impulsive noise, outliers, and inconsistent data.

This model explicitly allows considering the power fluctuation across range bins, especially for heavy-tailed clutter distributions.

The detection problem is written as a binary hypothesis test:

$$ \left\{ \begin{array}{ll} H_0: \mathbf{x}_i = \mathbf{z}_i & ; \quad \mathbf{x}_{i,k} = \mathbf{z}_{i,k} \quad k = 1 \dots K \\ H_1: \mathbf{x}_i = \alpha_i \mathbf{p}_i + \mathbf{z}_i & ; \quad \mathbf{x}_{i,k} = \mathbf{z}_{i,k} \quad k = 1 \dots K. \end{array} \right. \quad (4) $$

In (4) it is assumed that $K \ge m$ independent and identically distributed (i.i.d.) signal-free secondary data $\mathbf{x}_{i,k} \in \mathbb{C}^m$ are available under both hypotheses for background parameters estimation. We recall that $\mathbf{z}_{i,k} \stackrel{d}{=} \sqrt{\tau_{i,k}} \mathbf{c}_{i,k}$, $\mathbf{c}_{i,k} \sim \mathcal{CN}(\mathbf{0}, \mathbf{M}_{ii})$.

Conditionally to the unknown deterministic texture, the densities of $\mathbf{x}_i$ under $H_0$ and $H_1$ are Gaussian:

$$ \begin{aligned} p_{\mathbf{x}_i}(\mathbf{x}_i; H_0) &= \frac{1}{\pi^m \tilde{\sigma}_i^{2m} |\mathbf{M}_{ii}|} \exp \left( -\frac{\mathbf{x}_i^H \mathbf{M}_{ii}^{-1} \mathbf{x}_i}{\tilde{\sigma}_i^2} \right), \\ p_{\mathbf{x}_i}(\mathbf{x}_i; H_1) &= \frac{1}{\pi^m \tilde{\sigma}_i^{2m} |\mathbf{M}_{ii}|} \exp \left( -\frac{(\mathbf{x}_i - \alpha_i \mathbf{p}_i)^H \mathbf{M}_{ii}^{-1} (\mathbf{x}_i - \alpha_i \mathbf{p}_i)}{\tilde{\sigma}_i^2} \right), \end{aligned} \quad (5) $$

where $\tilde{\sigma}_i = \sigma_i \sqrt{\tau_i}$.

#### B. Data model for the two arrays

If we consider the two antennas, this model can be written more appropriately:

$$ \left\{ \begin{array}{ll} H_0: \mathbf{x} = \mathbf{z} & ; \quad \mathbf{x}_k = \mathbf{z}_k \quad k = 1 \dots K \\ H_1: \mathbf{x} = \mathbf{P} \alpha + \mathbf{z} & ; \quad \mathbf{x}_k = \mathbf{z}_k \quad k = 1 \dots K \end{array} \right. \quad (6) $$

where $\mathbf{x} = [\mathbf{x}_1 \ \mathbf{x}_2] \in \mathbb{C}^{2m}$ is the concatenation of the two received signals, $\alpha = \begin{bmatrix} \alpha_1 \\ \alpha_2 \end{bmatrix} \in \mathbb{C}^2$ is the vector of the target amplitudes and $\mathbf{z} = \begin{bmatrix} \mathbf{z}_1 \\ \mathbf{z}_2 \end{bmatrix} \in \mathbb{C}^{2m}$ is the additive clutter. The matrix $\mathbf{P} = \begin{bmatrix} \mathbf{p}_1 & \mathbf{0} \\ \mathbf{0} & \mathbf{p}_2 \end{bmatrix} \in \mathbb{C}^{2m \times 2}$ contains the steering vectors. $\{\mathbf{x}_k = \begin{bmatrix} \mathbf{x}_{1,k} \\ \mathbf{x}_{2,k} \end{bmatrix}\}_{k \in [1,K]} \in \mathbb{C}^{2m}$ for $K \ge 2m$ are i.i.d. signal-free secondary data.

This formulation allows considering the correlation between sensors of arrays. The covariance is $\mathbf{M} = \begin{bmatrix} \mathbf{M}_{11} & \mathbf{M}_{12} \\ \mathbf{M}_{21} & \mathbf{M}_{22} \end{bmatrix}$, with $\mathbf{M}_{ii}$ the block-covariance matrix of array $i$, and $\mathbf{M}_{ij} = \mathbf{M}_{ji}^H$ the cross-correlation block of array $i$ and $j$.

We further assume that the covariance is known, or estimated, up to two scalars, possibly different on each array. These scalars are conditioning the $\mathbf{M}_{ii}$ block-covariance, but also all the cross-correlations blocks associated with the array $i$.

We can therefore write:

$$ \tilde{\mathbf{C}} = \tilde{\mathbf{\Sigma}} \mathbf{M} \tilde{\mathbf{\Sigma}}, \quad (7) $$

with $\tilde{\mathbf{\Sigma}} = \begin{bmatrix} \tilde{\sigma}_1 \mathbf{I}_m & \mathbf{0} \\ \mathbf{0} & \tilde{\sigma}_2 \mathbf{I}_m \end{bmatrix} = \begin{bmatrix} \tilde{\sigma}_1 & 0 \\ 0 & \tilde{\sigma}_2 \end{bmatrix} \otimes \mathbf{I}_m$ the unknown diagonal matrix of scalars $\sigma_i$ and $\tau_i$. We remind that the
$\sigma_i$ parameter drives the partial homogeneity of the data (i.e. the difference in scale factor between the covariance matrices of the primary and secondary data) and $\tau_i$ drives the non-Gaussianity of the data (i.e. the power variation of the observations over time).

In this model, each array, although correlated, has a possibly distinct texture and an unknown scaling factor on the covariance matrix which may also be dissimilar. It would therefore become entirely possible to model, as an example, a first array whose observations would be Gaussian ($\tau_1 = 1$) and homogeneous ($\sigma_1 = 1$), and a second array whose observations would be K-distributed (for which $\tau_2$ is a realization of a Gamma distribution) and non-homogeneous ($\sigma_2 \neq 1$).

The PDFs under each hypothesis can be rewritten as:

$$ 
\begin{aligned}
p_{\mathbf{x}}(\mathbf{x}; H_0) &= \frac{1}{\pi^{2m} |\tilde{\mathbf{C}}|} \exp \left( -\mathbf{x}^H \tilde{\mathbf{C}}^{-1} \mathbf{x} \right), \\
p_{\mathbf{x}}(\mathbf{x}; H_1) &= \frac{1}{\pi^{2m} |\tilde{\mathbf{C}}|} \exp \left( -(\mathbf{x} - \mathbf{P}\boldsymbol{\alpha})^H \tilde{\mathbf{C}}^{-1} (\mathbf{x} - \mathbf{P}\boldsymbol{\alpha}) \right).
\end{aligned} 
\quad (8) $$

## IV. ROBUST DETECTORS

We discuss the derivation of detectors using the GLRT and Rao procedures. Following a two-step approach, the covariance matrix $\mathbf{M}$ will first be assumed known, and then replaced by an appropriate estimator.

### A. Detectors' derivation with $\mathbf{M}$ known (step-1)

1) Generalized Likelihood Ratio Test: The Generalized Likelihood Ratio Test (GLRT) design methodology proposes to solve the detection problem from the ratio of the probability densities function under $H_1$ and $H_0$, substituting the unknown parameters with their maximum likelihood estimates:

$$ L_G(\mathbf{x}) = \frac{\max_{\boldsymbol{\alpha}} \max_{\tilde{\boldsymbol{\Sigma}}} p_{\mathbf{x}}\left(\mathbf{x}; \boldsymbol{\alpha}, \tilde{\boldsymbol{\Sigma}}, H_1\right)}{\max_{\tilde{\boldsymbol{\Sigma}}} p_{\mathbf{x}}\left(\mathbf{x}; \tilde{\boldsymbol{\Sigma}}, H_0\right)}. \quad (9) $$

**Proposition IV.1.** The GLRT for the hypothesis test (6) is given by:

$$ L_G(\mathbf{x}) = \frac{\hat{\sigma}_{10} \hat{\sigma}_{20}}{\hat{\sigma}_{11} \hat{\sigma}_{21}} \quad (10) $$

where

$$ (\hat{\sigma}_{10} \hat{\sigma}_{20})^2 = \frac{\mathrm{Re}(\mathbf{x}_1^H \mathbf{M}_{12}^{-1} \mathbf{x}_2)}{m} + \sqrt{\frac{\mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1}{m} \frac{\mathbf{x}_2^H \mathbf{M}_{22}^{-1} \mathbf{x}_2}{m}}, $$

and

$$ (\hat{\sigma}_{11} \hat{\sigma}_{21})^2 = \frac{\mathrm{Re}[\mathbf{x}_1^H (\mathbf{M}_{12}^{-1} - \mathbf{D}_{12}^{-1}) \mathbf{x}_2]}{m} + \sqrt{\frac{\mathbf{x}_1^H (\mathbf{M}_{11}^{-1} - \mathbf{D}_{11}^{-1}) \mathbf{x}_1}{m} \frac{\mathbf{x}_2^H (\mathbf{M}_{22}^{-1} - \mathbf{D}_{22}^{-1}) \mathbf{x}_2}{m}} $$

$$ \text{with } \mathbf{D}^{-1} = \mathbf{M}^{-1} \mathbf{P} (\mathbf{P}^H \mathbf{M}^{-1} \mathbf{P})^{-1} \mathbf{P}^H \mathbf{M}^{-1}. $$

*Proof.* See Appendix A for a step-by-step derivation and Appendix B for some interesting equivalences. □

As $\mathbf{x}_i = \sqrt{\tau_i} \mathbf{c}_i$ under $H_0$, it is easily shown that this detector is texture independent (i.e. it has the texture-CFAR property).

This detection test will be called M-NMF-G in the following.

2) Rao test: The Rao test is obtained by exploiting the asymptotic efficiency of the ML estimate and expanding the likelihood ratio in the neighborhood of the estimated parameters [47]. A traditional approach for complex unknown parameters is to form a corresponding real-valued parameter vector and then use the real Rao test [48], [49]. Specifically, rewriting the detection problem (6) as:

$$ \begin{cases} H_0: \xi_R = 0, \xi_S \\ H_1: \xi_R \neq 0, \xi_S \end{cases} \quad (11) $$

where $\xi_R = [\mathrm{Re}(\boldsymbol{\alpha})^T \mathrm{Im}(\boldsymbol{\alpha})^T]^T$ is a $4 \times 1$ parameter vector and $\xi_S = [\tilde{\sigma}_1 \tilde{\sigma}_2]^T$ is a $2 \times 1$ vector of nuisance parameters, the Rao test for the problem (11) is:

$$ L_R(\mathbf{x}) = \frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \boldsymbol{\xi}_R, \boldsymbol{\xi}_S)}{\partial \boldsymbol{\xi}_R} \Bigg|_{\substack{\boldsymbol{\xi}_R = \hat{\boldsymbol{\xi}}_{R0} \\ \boldsymbol{\xi}_S = \hat{\boldsymbol{\xi}}_{S0}}} \left[ \mathbf{I}^{-1}(\hat{\boldsymbol{\xi}}_{R0}, \hat{\boldsymbol{\xi}}_{S0}) \right]_{\boldsymbol{\xi}_R \boldsymbol{\xi}_R} \frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \boldsymbol{\xi}_R, \boldsymbol{\xi}_S)}{\partial \boldsymbol{\xi}_R} \Bigg|_{\substack{\boldsymbol{\xi}_R = \hat{\boldsymbol{\xi}}_{R0} \\ \boldsymbol{\xi}_S = \hat{\boldsymbol{\xi}}_{S0}}} \quad (12) $$

The PDF $p_{\mathbf{x}}(\mathbf{x}; \boldsymbol{\xi}_R, \boldsymbol{\xi}_S)$ is given in (8) and parametrized by $\boldsymbol{\xi}_R$ and $\boldsymbol{\xi}_S$. $\hat{\boldsymbol{\xi}}_{R0}$ and $\hat{\boldsymbol{\xi}}_{S0}$ are the ML estimates of $\boldsymbol{\xi}_R$ and $\boldsymbol{\xi}_S$ under $H_0$. $\mathbf{I}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S)$ is the Fisher Information Matrix that can be partitioned as:

$$ \mathbf{I}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) = \begin{bmatrix} \mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_R}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) & \mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) \\ \mathbf{I}_{\boldsymbol{\xi}_S \boldsymbol{\xi}_R}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) & \mathbf{I}_{\boldsymbol{\xi}_S \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) \end{bmatrix}, $$

and we have:

$$ [\mathbf{I}^{-1}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S)]_{\boldsymbol{\xi}_R \boldsymbol{\xi}_R} = (\mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_R}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) - \mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) \mathbf{I}_{\boldsymbol{\xi}_S \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S))^{-1} \quad (13) $$

The following proposition can be finally stated:

**Proposition IV.2.** The Rao test is given by:

$$ L_R(\mathbf{x}) = 2 \mathbf{x}^H \hat{\mathbf{C}}_0^{-1} \mathbf{P} (\mathbf{P}^H \hat{\mathbf{C}}_0^{-1} \mathbf{P})^{-1} \mathbf{P}^H \hat{\mathbf{C}}_0^{-1} \mathbf{x} \quad (14) $$

where $\hat{\mathbf{C}}_0 = \hat{\boldsymbol{\Sigma}}_0 \mathbf{M} \hat{\boldsymbol{\Sigma}}_0$ and $\hat{\boldsymbol{\Sigma}}_0$ is defined in Appendix A.

*Proof.* See Appendix C. □

As for the GLRT, (14) is texture-CFAR. This detector will be referred to as M-NMF-R in the sequel.
## B. Covariance estimation and adaptive detectors (step-2)

In practice, the noise covariance matrix is unknown and estimated using the K available i.i.d. signal-free secondary data.

In Gaussian environment, in which PDFs are given by (8) with $\sigma_i = 1$ and $\tau_i = 1$, the MLE of M is the well-known Sample Covariance Matrix (SCM):

$$ \hat{\mathbf{M}}_{\text{SCM}} = \frac{1}{K} \sum_{k=1}^{K} \mathbf{x}_k \mathbf{x}_k^H \quad (15) $$

which is an unbiased and minimum variance estimator.

In the presence of outliers or a heavy-tailed distribution (as modeled by a mixture of scaled Gaussian), this estimator is no longer optimal or robust. This leads to a strong performance degradation.

**Proposition IV.3.** In MSG environment the MLE of M is given by:

$$ \hat{\mathbf{M}}_{2TYL} = \frac{1}{K} \sum_{k=1}^{K} \hat{\mathbf{T}}_{k}^{-1} \mathbf{x}_{k} \mathbf{x}_{k}^{H} \hat{\mathbf{T}}_{k}^{-1}, \quad (16) $$

where $\hat{\mathbf{T}}_k = \begin{bmatrix} \sqrt{\hat{\tau}_{1k}} & 0 \\ 0 & \sqrt{\hat{\tau}_{2k}} \end{bmatrix} \otimes \mathbf{I}_m$, $\hat{\tau}_{1k} = t_1 + \sqrt{\frac{t_1}{t_2}} t_{12}$, $\hat{\tau}_{2k} = t_2 + \sqrt{\frac{t_2}{t_1}} t_{12}$ and $t_1 = \frac{\mathbf{x}_{1,k}^H \hat{\mathbf{M}}_{11}^{-1} \mathbf{x}_{1,k}}{m}$, $t_2 = \frac{\mathbf{x}_{2,k}^H \hat{\mathbf{M}}_{22}^{-1} \mathbf{x}_{2,k}}{m}$, $t_{12} = \frac{\text{Re}(\mathbf{x}_{1,k}^H \hat{\mathbf{M}}_{12}^{-1} \mathbf{x}_{2,k})}{m}$.

*Proof.* The demonstration is provided in Appendix D. □

A key point is that this estimator is independent of the textures, i.e. the power variations, on each array. It can be thought of as a *multi-texture* generalization of [38].

From a practical point of view $\hat{\mathbf{M}}_{2TYL}$ is the solution of the recursive algorithm:

$$ \hat{\mathbf{T}}_{k}^{(n)} = \begin{bmatrix} \sqrt{\hat{\tau}_{1k}^{(n)}} & 0 \\ 0 & \sqrt{\hat{\tau}_{2k}^{(n)}} \end{bmatrix} \otimes \mathbf{I}_m, \quad (17) $$

$$ \hat{\mathbf{M}}_{2TYL}^{(n)} = \frac{1}{K} \sum_{k=1}^{K} \left( \hat{\mathbf{T}}_{k}^{(n-1)} \right)^{-1} \mathbf{x}_{k} \mathbf{x}_{k}^{H} \left( \hat{\mathbf{T}}_{k}^{(n-1)} \right)^{-1}, \quad (18) $$

where $n \in \mathbb{N}$ is the iteration number, whatever the initialization $\hat{\mathbf{M}}_{2TYL}^{(0)}$.

The convergence of recursive equations (17) and (18) in the estimation of $\hat{\mathbf{M}}_{2TYL}$ is illustrated in Figure 4 for 500 iterations and $\hat{\mathbf{M}}_{2TYL}^{(0)} = \mathbf{I}_{2m}$.

The relative difference between estimates decreases with the number of iterations. From iteration 60, the accuracy becomes limited by the simulation environment. In practice, it is not necessary to go to this limit, and we notice that from iteration 20 the relative deviation becomes lower than $10^{-6}$.

At last, the adaptive versions of the tests (10) and (14) will be simply obtained by replacing the covariance matrix M by an appropriate estimate: (15) or (16) according to

<figure><img src="image_5.png" /><figcaption>Fig. 4. Relative deviations of the estimates as a function of the number of iterations ($\hat{\mathbf{M}}_{2TYL}^{(0)} = \mathbf{I}_{2m}$).</figcaption></figure>

the environment. Those detectors will be referred as M-ANMF-$G_{SCM}$, M-ANMF-$R_{SCM}$, M-ANMF-$G_{TYL}$, and M-ANMF-$R_{TYL}$.

## V. NUMERICAL RESULTS ON SIMULATED DATA

### A. Performance assessment

Two correlation coefficients $\rho_1$ and $\rho_2$ ($0 < \rho_1, \rho_2 < 1$), are used in the construction of a *speckle* covariance matrix model defined as:

$$ [\mathbf{M}]_{jl} = \beta \rho_1^{|j_1-l_1|} \times \rho_2^{|j_2-l_2|}, $$

where $j, l \in [1, 2m]$ are sensor numbers of coordinates ($j_1, j_2$) and ($l_1, l_2$) respectively and $\beta$ is a scale factor.

Thus, denoting $\mathbf{M} = \begin{bmatrix} \mathbf{M}_{11} & \mathbf{M}_{12} \\ \mathbf{M}_{21} & \mathbf{M}_{22} \end{bmatrix}$, $\mathbf{M}_{11}$ and $\mathbf{M}_{22}$ are the covariance matrices of array 1 and 2 and $\mathbf{M}_{12} = \mathbf{M}_{21}^H$ is a cross-correlation block. In the FLS context, the block $\mathbf{M}_{11}$ is weakly correlated and Toeplitz structured (close to the identity matrix). $\mathbf{M}_{22}$ is also Toeplitz but more strongly correlated (due to the wider transmission beam). The cross-correlation blocks could be considered null under the uncorrelated arrays assumption. In our case, these will be different from zero because the arrays cross each other in their central part. This results in the general structure displayed in Figure 5, where we visually show the adequacy of this model with the SCM covariance estimator established on real data.

We choose $\sigma_i = 1$. $\tau_i = 1$ for Gaussian clutter and $\tau_i \sim \text{Gam}(\nu, 1/\nu)$ with $\nu = 0.5$ for impulsive non-Gaussian K-distributed data (that is, the texture variables follow a gamma distribution with shape parameter 0.5 and scale parameter 2). The PFA-threshold curves are established on the basis of 1000 realizations of random vectors. The detection probabilities are statistically estimated by adding a synthetic narrow band far field point target with identical amplitudes on the arrays ($\alpha_1 = \alpha_2$). 10000 Monte Carlo iterations are performed and 150 target amplitudes are evaluated.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Sonar arrays (left): two uniform linear arrays intersect in their centers. Empirical covariance matrix on real sonar data shown in dB (center): The SCM is estimated on a homogeneous area between 265 m and 328 m. Cross-correlation blocks are non-zero. The covariance matrix model used for the data simulation is shown in dB (right): $\beta = 3 \times 10^{-4}$, $\rho_1 = 0.4$, $\rho_2 = 0.9$.</figcaption></figure>

## B. Benchmark tests

The GLRT for a single array in a partially homogeneous Gaussian environment when $M_{ii}$ is known is the Normalized Matched Filter [11]:

$$ \text{NMF}_i(\mathbf{x}_i) = \frac{|\mathbf{p}_i^H \mathbf{M}_{ii}^{-1} \mathbf{x}_i|^2}{(\mathbf{p}_i^H \mathbf{M}_{ii}^{-1} \mathbf{p}_i) (\mathbf{x}_i^H \mathbf{M}_{ii}^{-1} \mathbf{x}_i)}. \quad (19) $$

Adaptive versions are obtained by substituting the covariance matrix by a suitable estimate [35] and will be referred to as ANMF$_{SCM}$i in Gaussian case, or ANMF$_{TYL}$i for the non-Gaussian case.

When the two arrays are considered, in the very favorable case of a Gaussian homogeneous environment where the covariance matrix C is perfectly known, the GLRT is:

$$ \text{MIMO-MF}(\mathbf{x}) = \mathbf{x}^H \mathbf{C}^{-1} \mathbf{P} (\mathbf{P}^H \mathbf{C}^{-1} \mathbf{P})^{-1} \mathbf{P}^H \mathbf{C}^{-1} \mathbf{x}. \quad (20) $$

This is the MIMO Optimum Gaussian Detector (R-MIMO OGD) in [50], which is a multi-array generalization of the Matched Filter test. One can note the very strong similarity with (14). Its adaptive version is MIMO-AMF$_{SCM}$. It seems useful to specify that this detector is relevant only in a Gaussian and perfectly homogeneous environment. Especially, exactly as in the single-array case [31], the covariance estimator (16) is defined up to a constant scalar. Detectors (10) and (14) are invariant when the covariance is changed by a scale factor. This is a direct result of the partial homogeneity assumption. This is not the case for (20), thus the adaptive MIMO-AMF$_{TYL}$ version is not relevant and will not be used in the following.

These tests will be considered as performance bounds for the proposed detectors: M-NMF-G and M-NMF-R when M is known, M-ANMF-G$_{SCM}$, M-ANMF-R$_{SCM}$, M-ANMF-G$_{TYL}$ and M-ANMF-R$_{TYL}$ otherwise.

## C. Performance in Gaussian clutter

We have shown that all developed detectors are texture-CFAR. Unfortunately, the matrix CFAR property (distribution of the test keeps identical whatever the true covariance matrix) is much more difficult to show. Therefore, we propose to perform a study on simulated data to check this matrix CFAR property. Figure 6 experimentally demonstrates the CFAR behavior of the detectors with respect to the covariance matrix in a Gaussian environment.

On the left side we represent the false alarm probability of Rao's detector (for known and estimated M) as a function of the detection threshold. On the right, these curves are plotted for the GLRT detector.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Pfa-threshold curves of the M-(A)NMF-R$_{SCM}$ and M-(A)NMF-G$_{SCM}$ detectors in Gaussian environment. Left: the Rao detector for known M (blue), for M estimated from $2 \times 2m$ secondary data (red), the OGD detector or Matched Filter (black). Right: The GLRT detector for known M (blue), and estimated based on $2 \times 2m$ secondary data (red).</figcaption></figure>

The deviation of the adaptive detectors comes from the covariance estimation process: the red curve converges to the blue one when the number of secondary data used for the estimation of M increases.

The markers represent curves for different covariances. In solid line $\beta = 3 \times 10^{-4}$, $\rho_1 = 0.4$, $\rho_2 = 0.9$. The cross markers are established from $\beta = 1$, $\rho_1 = 0.95$, $\rho_2 = 0.95$. Circular markers are for $\beta = 100$, $\rho_1 = 0.1$, $\rho_2 = 0.1$ and null anti-diagonal blocks. For very distinct covariances, the superposition of curves and markers underlines the matrix-CFAR property of the detectors.

Figure 7 illustrates the value of merging arrays with a cross-shaped geometry. The detection probabilities of NMF 1, NMF 2, and M-NMF-R are plotted as a function of azimuth ($\theta_1$) and elevation ($\theta_2$) angles at fixed SNR.

In the first angular dimension, which depends on the array considered, each ULA has a zone of least detectability close to 0°. This is due to the correlated nature of clutter. The decrease in detection probabilities from +/-60° to 0° is a function of the correlation level ($\rho_1$ or $\rho_2$). In the second dimension, the detection performances are invariant. Thus, the poor detection performance near 0° propagates in a whole direction of space: "vertically" for NMF 1 and "horizontally" for NMF 2. As an illustration, the probability of detection of NMF 1 in $\theta_1 = 0°$ whatever $\theta_2$ is $PD = 0.17$, and the probability of detection of NMF 2 in $\theta_2 = 0°$ for all $\theta_1$ is 0.03.

The M-NMF-R detector spatially minimizes this area of least detectability. In this case, the probability of detection at $\theta_1 = 0°$ becomes greater than 0.8 for $|\theta_2| > 10.5°$ and
<figure><img src="image_8.png" /><figcaption>Fig. 7. Probability of detection in Gaussian environment for known covariance matrix **M** as a function of $\theta_1$ and $\theta_2$ ($\beta = 3 \times 10^{-4}$, $\rho_1 = 0.4$, $\rho_2 = 0.9$, SNR = -12 dB and PFA = $10^{-2}$). Top left: NMF 1. Top right: NMF 2. Bottom left: M-NMF-R.</figcaption></figure>

greater than 0.5 for $\theta_2 = 0^\circ$ and $|\theta_1| > 24^\circ$. The *Rigde* clutter is minimized.

The probability of detection (PD) is plotted as a function of the signal-to-noise ratio (SNR) in Figure 8 when **M** is known.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Probability of detection in Gaussian environment for known $\mathbf{M}$ ($\beta = 3 \times 10^{-4}$, $\rho_1 = 0.4$, $\rho_2 = 0.9$ and PFA = $10^{-2}$).</figcaption></figure>

The MIMO-MF detector (shown in black), which requires a perfect knowledge of the covariance, is logically the most efficient. The single array NMF detectors (NMF 1 and 2, in purple and green) have comparable and the lowest performance. The M-NMF-I detector (in blue) between these curves assumes antenna independence. The proposed M-NMF-G (red) and M-NMF-R (yellow) detectors are both equivalent and superior to the M-NMF-I (0.2 dB at

$PD = 0.8$) and much more efficient than NMF 1 (2.5 dB) and NMF 2 tests (2 dB). The difference with MIMO-MF is slight, around 0.2 dB at $PD = 0.8$. Detection performances by the Rao approach are slightly superior to those obtained by the GLRT.

Figure 9 compares the performance of the tests in their adaptive versions. The MIMO-MF curve is shown in black dotted lines for illustrative purposes only (as it assumes the covariance known). It can be seen that when the SCM is evaluated on $2 \times 2m = 256$ secondary data, the covariance estimation leads to a loss of 3 dB between the performances of the MIMO-MF test and its adaptive version MIMO-AMF, as expected from the Reed-Mallett-Brennan's rule [51].

<figure><img src="image_10.png" /><figcaption>Fig. 9. Probability of detection in Gaussian environment for unknown $\mathbf{M}$ (PFA = $10^{-2}$, SCMs are estimated based on $2 \times 2m$ secondary data).</figcaption></figure>

In their adaptive versions, the proposed detectors offer equivalent performances to the MIMO-AMF detectors while offering additional robustness and flexibility conditions on a possible lack of knowledge of the covariance (estimated to be within two scale factors). The gain compared to single array detectors, although reduced compared to the case where **M** is known, remains favorable and of the order of 0.5 to 1 dB at $PD = 0.8$. In other words, for an SNR of -11 dB, $PD = 0.75$ for ANMF 1 and 0.85 for M-ANMF-G$_{SCM}$ (or M-ANMF-R$_{SCM}$).

### D. Performance in impulsive non-Gaussian clutter

PFA-threshold curves in *K*-distributed environment are displayed in Figure 10 in order to characterize the matrix-CFAR behavior.

The detectors based on the new covariance matrix estimator systematically have lower detection thresholds to those obtained with the SCM. While optimal for a Gaussian clutter, the MIMO-MF detector is no longer suitable. The marker overlays highlight the matrix-CFAR behavior of the Rao detector in a non-Gaussian environment. The
<figure><img src="image_11.png" /><figcaption>Fig. 10. Pfa-threshold curves in non-Gaussian K-distributed environment ($\nu = 0.5$). Left: the Rao detector for known M (blue), for M estimated with the SCM (red), M estimated with the new estimator (yellow), the OGD detector or *Matched Filter* (black). Right: The GLRT detector for known M (blue), M estimated based on the SCM (red), based on Tyler's extension (yellow).</figcaption></figure>

case of the GLRT detector is much less obvious: at this stage, it seems complicated to consider that the false alarm curves are strictly independent of the covariance matrix.

The detection performances are compared in Figure 11.

<figure><img src="image_12.png" /><figcaption>Fig. 11. Probability of detection in K-distributed environment for unknown M ($\nu = 0.5$, $PFA = 10^{-2}$, $2 \times 2m$ secondary data). In solid lines, the detectors are built with the SCM estimate. In dotted lines, the detectors are built with Tyler's (single array) or the new (dual arrays) estimator.</figcaption></figure>

For the Rao and GLRT detectors, using the estimator (16) in an impulsive non-Gaussian environment induces a gain in detection performance of the order of 1.2 dB at $PD = 0.8$ with the SCM. Compared to the best competitors the improvement is 3 dB at $PD = 0.8$.

# VI. EXPERIMENTAL DATA WITH REAL TARGET

The actual sonar measurements were collected during an acquisition campaign in May 2022 in La Ciotat bay, Mediterranean Sea. The experimental conditions are described in II-B, and the dataset of interest is shown in Figure 3. Echoes from the true target (right in Figure 2) are observed at the 2040th range cell.

## A. Illustrations of detection test outputs

Figure 12 provides examples of real data outputs from conventional ANMF detectors based on a single array. The real target is observed in range bin 2040 (143 m away from the sonar), at angular bin $\theta_1 = 26$ ($-12.4^\circ$ azimuth), and angular bin $\theta_2 = 37$ ($8.6^\circ$ elevation).

<figure><img src="image_13.png" /><figcaption>Fig. 12. ANMF$_{SCM}$ detector outputs on real sonar data. A target is observed on array 1 (left) at coordinates (26, 2040) and on array 2 (right) at coordinates (37, 2040). The SCM is built from 256 secondary data.</figcaption></figure>

Figure 13 shows the outputs of the Rao and GLRT detectors applied simultaneously to both arrays at the specific range bin 2040. These subfigures are not directly comparable with each other or with Figure 12. The target is perfectly located in azimuth and elevation.

<figure><img src="image_14.png" /><figcaption>Fig. 13. Outputs of the M-ANMF-R$_{SCM}$ (left) and M-ANMF-G$_{SCM}$ (right) detectors. The real target is at coordinates (37, 26). The SCM is built from 256 secondary data.</figcaption></figure>

## B. Robustness to training data corruption

Previously, we assumed the availability of a set of $K \ge 2m$ i.i.d. secondary data, supposed free of signal components and sharing statistic properties with the noise in the cell under test. In practice, such data can be obtained by processing samples in spatial proximity to the range bin being evaluated, and the absence of signals is not always verified or checkable. In particular, this assumption is no longer valid if another target is present in these secondary data.
Figure 14 illustrates the robustness of the covariance matrix estimators to data corruption. A synthetic target is added 100 samples away from the real target and contaminates the secondary dataset.

<figure><img src="image_15.png" /><figcaption>Fig. 14. Rao detector outputs based on the SCM (left) or Tyler (right) with secondary data corruption. The covariance matrix estimators are based on 256 secondary data.</figcaption></figure>

On the left side, the output of M-ANMF-R<sub>SCM</sub> is strongly degraded. The SCM is not robust to the presence of outliers. The target is hardly discernible, and the maximum is wrongly located.

Under the same conditions, the response of the M-ANMF-R<sub>TYL</sub> detector is visualized in the right part. Although degraded compared to Figure 13, the behavior is still largely usable. The target presence is easily identified, and the new estimator (16) is much more resistant to data contamination.

## VII. CONCLUSIONS

In this paper, we considered the problem of adaptive point target detection by a correlated multi-arrays Mills Cross sonar system. Using a 2-step approach, we first derived two new detectors that are robust to differences in target amplitudes and to unknown scaling factors on the covariances. Subsequently, we have introduced an innovative estimator of the covariance matrix suitable to any non-Gaussian MSG environment. By these very general assumptions, the framework of the study can therefore concern systems with co-located or remote arrays.

Experimental results show that the detection performance is up to 3 dB better than conventional approaches. The detectors cover a larger detection area and are particularly robust to spikes, impulsive noise, and data contamination.

Future work will focus on establishing a theoretical demonstration of the matrix-CFAR behavior of these detectors, and on generalizing solutions for different numbers and geometries of arrays.

## ACKNOWLEDGMENTS

This work has been done thanks to the facilities offered by the Univ. Savoie Mont Blanc - CNRS/IN2P3 MUST computing center.

## APPENDIX A
GLRT'S DERIVATION

For the following, and for ease of reading, the punctuation mark “tilde” will be omitted. Thus $\Sigma$ and $\hat{\sigma}_i$ will

simply be denoted as $\Sigma$ and $\sigma_i$. The same is true for their respective estimates.

### A. Maximum Likelihood Estimator of $\Sigma$ under $H_0$

The MLE $\hat{\Sigma}_0$ is derived from the log-likelihood function:

$$ \ln p_{\mathbf{x}}(\mathbf{x}; \mathbf{\Sigma}, H_0) = -mL \ln \pi + 2 \ln |\mathbf{\Sigma}^{-1}| - \ln |\mathbf{M}| - (\mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x}) $$

whose derivative relative to $\mathbf{\Sigma}^{-1}$ is:

$$ \frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \mathbf{\Sigma}, H_0)}{\partial \mathbf{\Sigma}^{-1}} = 2 \frac{\partial \ln |\mathbf{\Sigma}^{-1}|}{\partial \mathbf{\Sigma}^{-1}} - \frac{\partial \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x}}{\partial \mathbf{\Sigma}^{-1}} $$

Knowing that $([52] (57)) \frac{\partial \ln |\mathbf{\Sigma}^{-1}|}{\partial \mathbf{\Sigma}^{-1}} = \mathbf{\Sigma}$, and $([52] (82))$

$$ \frac{\partial \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x}}{\partial \mathbf{\Sigma}^{-1}} = 2 \text{Re} (\mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x} \mathbf{x}^H) \text{ we have:} $$

$$ \frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \mathbf{\Sigma}, H_0)}{\partial \mathbf{\Sigma}^{-1}} = 2\mathbf{\Sigma} - 2\text{Re}(\mathbf{M}^{-1}\mathbf{\Sigma}^{-1}\mathbf{x}\mathbf{x}^H), $$

which leads to:

$$ \hat{\Sigma}_0 = \text{Re} (\mathbf{M}^{-1} \hat{\Sigma}_0^{-1} \mathbf{x} \mathbf{x}^H). \quad (21) $$

Expanding this matrix product with $\hat{\Sigma}_0 = [\hat{\sigma}_{10} \mathbf{I}_m \quad \mathbf{0} \quad \hat{\sigma}_{20} \mathbf{I}_m]$, we have:

$$ \hat{\sigma}_{10} \mathbf{I}_m = \text{Re} \left( \mathbf{M}_{11}^{-1} \frac{\mathbf{x}_1 \mathbf{x}_1^H}{\hat{\sigma}_{10}} + \mathbf{M}_{12}^{-1} \frac{\mathbf{x}_2 \mathbf{x}_1^H}{\hat{\sigma}_{20}} \right), $$

and using the trace operator:

$$ m \hat{\sigma}_{10} = \text{tr} \left[ \text{Re} \left( \mathbf{M}_{11}^{-1} \frac{\mathbf{x}_1 \mathbf{x}_1^H}{\hat{\sigma}_{10}} + \mathbf{M}_{12}^{-1} \frac{\mathbf{x}_2 \mathbf{x}_1^H}{\hat{\sigma}_{20}} \right) \right] \\ = \frac{\mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1}{\hat{\sigma}_{10}} + \text{Re} \left( \frac{\mathbf{x}_1^H \mathbf{M}_{12}^{-1} \mathbf{x}_2}{\hat{\sigma}_{20}} \right). $$

We then have:

$$ \hat{\sigma}_{10} = \frac{1}{\hat{\sigma}_{10}} \frac{\mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1}{m} + \frac{1}{\hat{\sigma}_{20}} \frac{\text{Re}(\mathbf{x}_1^H \mathbf{M}_{12}^{-1} \mathbf{x}_2)}{m}, \quad (22) $$

and

$$ \hat{\sigma}_{20} = \frac{1}{\hat{\sigma}_{10}} \frac{\text{Re}(\mathbf{x}_2^H \mathbf{M}_{21}^{-1} \mathbf{x}_1)}{m} + \frac{1}{\hat{\sigma}_{20}} \frac{\mathbf{x}_2^H \mathbf{M}_{22}^{-1} \mathbf{x}_2}{m}. \quad (23) $$

Denoting $a_1 = \frac{\mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1}{m}$, $a_{12} = \frac{\text{Re}(\mathbf{x}_1^H \mathbf{M}_{12}^{-1} \mathbf{x}_2)}{m}$, $a_2 = \frac{\mathbf{x}_2^H \mathbf{M}_{22}^{-1} \mathbf{x}_2}{m}$ and using (22) and (23):

$$ \left\{ \begin{array}{l} \hat{\sigma}_{10}^2 = a_1 + \frac{\hat{\sigma}_{10}}{\hat{\sigma}_{20}} a_{12} \\ \hat{\sigma}_{20}^2 = \frac{\hat{\sigma}_{20}}{\hat{\sigma}_{10}} a_{12} + a_2, \end{array} \right. \quad (24) $$

or:

$$ \left\{ \begin{array}{l} \hat{\sigma}_{10}^2 = a_1 + \frac{\hat{\sigma}_{10}}{\hat{\sigma}_{20}} a_{12} \\ \hat{\sigma}_{10}^2 = \frac{\hat{\sigma}_{10}}{\hat{\sigma}_{20}} a_{12} + \frac{\hat{\sigma}_{10}^2}{\hat{\sigma}_{20}^2} a_2. \end{array} \right. $$
By equalization of right-hand terms:

$$a_1 = \frac{\hat{\sigma}_{10}^2}{\hat{\sigma}_{20}^2} a_2,$$

and keeping the positive solution: $\frac{\hat{\sigma}_{10}}{\hat{\sigma}_{20}} = \sqrt{\frac{a_1}{a_2}}$, we obtain from (24):

$$\left\{ \begin{array}{l} \hat{\sigma}_{10}^2 = a_1 + \sqrt{\frac{a_1}{a_2}} a_{12} \\ \\ \hat{\sigma}_{20}^2 = a_2 + \sqrt{\frac{a_2}{a_1}} a_{12}, \end{array} \right. \qquad (25)$$

and

$$\hat{\sigma}_{10}^2 \hat{\sigma}_{20}^2 = (\sqrt{a_1} a_2 + a_{12})^2. \qquad (26)$$

## B. Maximum Likelihood Estimator of $\alpha$ under $H_1$

The ML estimates $\hat{\alpha} = [\hat{\alpha}_1 \hat{\alpha}_2]^T$ is found minimizing $(x - P\alpha)^H C^{-1} (x - P\alpha)$ with respect to $\alpha$ as ([53] (15.50)):

$$\begin{aligned} \hat{\alpha} &= (P^H C^{-1} P)^{-1} P^H C^{-1} x, \\ &= (P^H \Sigma^{-1} M^{-1} \Sigma^{-1} P)^{-1} P^H \Sigma^{-1} M^{-1} \Sigma^{-1} x. \end{aligned} \qquad (27)$$

## C. Maximum Likelihood Estimator of $\Sigma$ under $H_1$

The derivative of the log-likelihood function under $H_1$ with respect to $\Sigma^{-1}$ is:

$$\frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \hat{\alpha}, \mathbf{\Sigma}, H_1)}{\partial \mathbf{\Sigma}^{-1}} = 2 \frac{\partial \ln |\mathbf{\Sigma}^{-1}|}{\partial \mathbf{\Sigma}^{-1}} - \frac{\partial (\mathbf{x} - \mathbf{P}\hat{\alpha})^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} (\mathbf{x} - \mathbf{P}\hat{\alpha})}{\partial \mathbf{\Sigma}^{-1}}.$$

Furthermore:

$$\begin{aligned} (\mathbf{x} - \mathbf{P}\hat{\alpha})^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} (\mathbf{x} - \mathbf{P}\hat{\alpha}) &= \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x} - \hat{\alpha}^H \mathbf{P}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x} \\ &= \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x} - \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{P} \\ &\quad (\mathbf{P}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{P})^{-1} \mathbf{P}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x}. \end{aligned}$$

If $\mathbf{Z} = \begin{bmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{bmatrix}$, we then notice that $\mathbf{\Sigma}^{-1}\mathbf{P} = \mathbf{P}\mathbf{Z}^{-1}$ and $\mathbf{P}^H\mathbf{\Sigma}^{-1} = \mathbf{Z}^{-1}\mathbf{P}^H$. Thus:

$$\begin{aligned} & \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{P} (\mathbf{P}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{P})^{-1} \\ & \qquad \mathbf{P}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x} \\ & = \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{P}\mathbf{Z}^{-1} (\mathbf{Z}^{-1}\mathbf{P}^H \mathbf{M}^{-1}\mathbf{P}\mathbf{Z}^{-1})^{-1} \\ & \qquad \mathbf{Z}^{-1}\mathbf{P}^H \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x} \\ & = \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{P} (\mathbf{P}^H \mathbf{M}^{-1} \mathbf{P})^{-1} \mathbf{P}^H \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x}. \end{aligned}$$

By denoting $\mathbf{D}^{-1} = \mathbf{M}^{-1}\mathbf{P}(\mathbf{P}^H\mathbf{M}^{-1}\mathbf{P})^{-1}\mathbf{P}^H\mathbf{M}^{-1}$:

$$\begin{aligned} (\mathbf{x} - \mathbf{P}\hat{\alpha})^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} (\mathbf{x} - \mathbf{P}\hat{\alpha}), \\ = \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x} - \mathbf{x}^H \mathbf{\Sigma}^{-1} \mathbf{D}^{-1} \mathbf{\Sigma}^{-1} \mathbf{x}. \end{aligned}$$

So:

$$\frac{\partial (\mathbf{x} - \mathbf{P}\hat{\alpha})^H \mathbf{\Sigma}^{-1} \mathbf{M}^{-1} \mathbf{\Sigma}^{-1} (\mathbf{x} - \mathbf{P}\hat{\alpha})}{\partial \mathbf{\Sigma}^{-1}} = 2\mathrm{Re} \left[ (\mathbf{M}^{-1} - \mathbf{D}^{-1}) \mathbf{\Sigma}^{-1} \mathbf{x} \mathbf{x}^H \right].$$

Finally:

$$\frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \hat{\alpha}, \mathbf{\Sigma}, H_1)}{\partial \mathbf{\Sigma}^{-1}} = 2\mathbf{\Sigma} - 2\mathrm{Re} \left[ (\mathbf{M}^{-1} - \mathbf{D}^{-1}) \mathbf{\Sigma}^{-1} \mathbf{x} \mathbf{x}^H \right].$$

The minimum is given by:

$$\hat{\mathbf{\Sigma}}_1 = \mathrm{Re} \left[ (\mathbf{M}^{-1} - \mathbf{D}^{-1}) \hat{\mathbf{\Sigma}}_1^{-1} \mathbf{x} \mathbf{x}^H \right]. \qquad (28)$$

The rest is identical to paragraph A-A.

## D. Expression of the GLRT

As $\mathbf{x}^H \hat{\mathbf{\Sigma}}_0^{-1} \mathbf{M}^{-1} \hat{\mathbf{\Sigma}}_0^{-1} \mathbf{x}$ is a real positive scalar, we have:

$$\begin{aligned} \mathbf{x}^H \hat{\mathbf{\Sigma}}_0^{-1} \mathbf{M}^{-1} \hat{\mathbf{\Sigma}}_0^{-1} \mathbf{x} &= \mathrm{Re} \left[ \mathrm{tr} \left( \hat{\mathbf{\Sigma}}_0^{-1} \mathbf{M}^{-1} \hat{\mathbf{\Sigma}}_0^{-1} \mathbf{x} \mathbf{x}^H \right) \right], \\ &= \mathrm{tr} \left[ \hat{\mathbf{\Sigma}}_0^{-1} \mathrm{Re} \left( \mathbf{M}^{-1} \hat{\mathbf{\Sigma}}_0^{-1} \mathbf{x} \mathbf{x}^H \right) \right], \\ &= mL. \end{aligned}$$

So the two PDF take the form:

$$p_{\mathbf{x}}(\mathbf{x}; \hat{\mathbf{\Sigma}}_0, H_0) = \frac{1}{\pi^{mL} |\hat{\mathbf{\Sigma}}_0|^2 |\mathbf{M}|} \exp(-mL),$$

$$p_{\mathbf{x}}(\mathbf{x}; \hat{\alpha}, \hat{\mathbf{\Sigma}}_1, H_1) = \frac{1}{\pi^{mL} |\hat{\mathbf{\Sigma}}_1|^2 |\mathbf{M}|} \exp(-mL).$$

The GLRT is expressed as:

$$L_G(\mathbf{x}) = \frac{\frac{1}{|\hat{\mathbf{\Sigma}}_1|^2}}{\frac{1}{|\hat{\mathbf{\Sigma}}_0|^2}} = \frac{|\hat{\mathbf{\Sigma}}_0|^2}{|\hat{\mathbf{\Sigma}}_1|^2} \qquad (29)$$

which can be thought of as a generalized variance ratio.

As $\hat{\mathbf{\Sigma}}_0$ and $\hat{\mathbf{\Sigma}}_1$ are diagonal matrices, $|\hat{\mathbf{\Sigma}}_0| = \prod_{i=1}^2 \hat{\sigma}_{i0}^m$ and $|\hat{\mathbf{\Sigma}}_1| = \prod_{i=1}^2 \hat{\sigma}_{i1}^m$ which leads to (10).

# APPENDIX B

## SPECIAL CASES OF THE GLRT DETECTOR

The detector (10) can be declined in already-known cases.

### A. Single array case

For a single array, we have $\hat{\sigma}_{10}^2 = \frac{\mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1}{m}$ and $\hat{\sigma}_{11}^2 = \frac{1}{m} \left( \mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1 - \frac{|\mathbf{p}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1|^2}{\mathbf{p}_1^H \mathbf{M}_{11}^{-1} \mathbf{p}_1} \right)$.
Replacing into the likelihood $L_G(\mathbf{x}) = \frac{|\hat{\Sigma}_0|^2}{|\hat{\Sigma}_1|^2} = \frac{\hat{\sigma}_{10}^{2m}}{\hat{\sigma}_{11}^{2m}}$ leads to:

$$L_G(\mathbf{x})^{1/m} = \frac{\hat{\sigma}_{10}^2}{\hat{\sigma}_{11}^2} = \frac{\frac{\mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1}{m}}{\frac{1}{m} \left( \mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1 - \frac{|\mathbf{p}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1|^2}{\mathbf{p}_1^H \mathbf{M}_{11}^{-1} \mathbf{p}_1} \right)} \\ = \left( 1 - \frac{|\mathbf{p}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1|^2}{(\mathbf{p}_1^H \mathbf{M}_{11}^{-1} \mathbf{p}_1) (\mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1)} \right)^{-1}.$$

Defining $l_G(\mathbf{x}) = \frac{L_G(\mathbf{x})^{1/m} - 1}{L_G(\mathbf{x})^{1/m}} = 1 - L_G(\mathbf{x})^{-1/m}$, we obtain the well-known NMF detector [11]:

$$l_G(\mathbf{x}) = \frac{|\mathbf{p}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1|^2}{(\mathbf{p}_1^H \mathbf{M}_{11}^{-1} \mathbf{p}_1) (\mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1)}. \quad (30)$$

### B. Uncorrelated arrays case

When the two arrays are fully uncorrelated, we have $\hat{\sigma}_{i0}^2 = \frac{\mathbf{x}_i^H \mathbf{M}_{ii}^{-1} \mathbf{x}_i}{m}$ and $\hat{\sigma}_{i1}^2 = \frac{1}{m} \left( \mathbf{x}_i^H \mathbf{M}_{ii}^{-1} \mathbf{x}_i - \frac{|\mathbf{p}_i^H \mathbf{M}_{ii}^{-1} \mathbf{x}_i|^2}{\mathbf{p}_i^H \mathbf{M}_{ii}^{-1} \mathbf{p}_i} \right)$.

We obtain:

$$L_G(\mathbf{x}) = \frac{|\hat{\Sigma}_0|^2}{|\hat{\Sigma}_1|^2} = \frac{\prod_{i=1}^{2} \hat{\sigma}_{i0}^{2m}}{\prod_{i=1}^{2} \hat{\sigma}_{i1}^{2m}} = \prod_{i=1}^{2} \left[ \frac{\hat{\sigma}_{i0}^2}{\hat{\sigma}_{i1}^2} \right]^m, \quad (31)$$

$$= \prod_{i=1}^{2} \left[ 1 - \frac{|\mathbf{p}_i^H \mathbf{M}_{ii}^{-1} \mathbf{x}_i|^2}{(\mathbf{p}_i^H \mathbf{M}_{ii}^{-1} \mathbf{p}_i) (\mathbf{x}_i^H \mathbf{M}_{ii}^{-1} \mathbf{x}_i)} \right]^{-m}. \quad (32)$$

This corresponds to the MIMO ANMF detector on independent arrays presented in [40].

### C. $\Sigma = \sigma I_{2m}$ case

When $\Sigma$ is the identity matrix up to a scalar factor, $\sigma_1 = \sigma_2$, whose estimators are renamed $\hat{\sigma}_0$ under $H_0$ and $\hat{\sigma}_1$ under $H_1$.

$$L_G(\mathbf{x}) = \frac{|\hat{\Sigma}_0|^2}{|\hat{\Sigma}_1|^2}, \quad \text{with} \quad \hat{\Sigma}_0 = \hat{\sigma}_0 I_{mL} \quad \text{and} \quad \hat{\Sigma}_1 = \hat{\sigma}_1 I_{mL} \\ \frac{\hat{\sigma}_0^{2mL}}{\hat{\sigma}_1^{2mL}}.$$

From $\hat{\Sigma}_0 = \text{Re}(\mathbf{M}^{-1}\hat{\Sigma}_0^{-1}\mathbf{x}\mathbf{x}^H)$, we have the following relations:

$$\text{tr}(\hat{\Sigma}_0) = \text{tr}[\text{Re}(\mathbf{M}^{-1}\hat{\Sigma}_0^{-1}\mathbf{x}\mathbf{x}^H)],$$

$$\hat{\sigma}_0 \text{tr}(I_{mL}) = \frac{1}{\hat{\sigma}_0} \text{Re}[\text{tr}(\mathbf{M}^{-1}\mathbf{x}\mathbf{x}^H)],$$

$$\hat{\sigma}_0^2 = \frac{\mathbf{x}_1^H \mathbf{M}_{11}^{-1} \mathbf{x}_1}{mL}, \text{ as } \mathbf{M}^{-1} \text{ is positive definite.}$$

Identically, we have:

$$\hat{\sigma}_1^2 = \frac{(\mathbf{x} - \mathbf{P}\alpha)^H \mathbf{M}^{-1} (\mathbf{x} - \mathbf{P}\alpha)}{mL}, \\ = \frac{\mathbf{x}^H \mathbf{M}^{-1} \mathbf{x} - \mathbf{x}^H \mathbf{M}^{-1} \mathbf{P} (\mathbf{P}^H \mathbf{M}^{-1} \mathbf{P})^{-1} \mathbf{P}^H \mathbf{M}^{-1} \mathbf{x}}{mL}.$$

$$L_G(\mathbf{x})^{1/mL} = \frac{\hat{\sigma}_0^2}{\hat{\sigma}_1^2}, \\ = \frac{\mathbf{x}^H \mathbf{M}^{-1} \mathbf{x}}{\mathbf{x}^H \mathbf{M}^{-1} \mathbf{x} - \mathbf{x}^H \mathbf{M}^{-1} \mathbf{P} (\mathbf{P}^H \mathbf{M}^{-1} \mathbf{P})^{-1} \mathbf{P}^H \mathbf{M}^{-1} \mathbf{x}} \\ = \left( 1 - \frac{\mathbf{x}^H \mathbf{M}^{-1} \mathbf{P} (\mathbf{P}^H \mathbf{M}^{-1} \mathbf{P})^{-1} \mathbf{P}^H \mathbf{M}^{-1} \mathbf{x}}{\mathbf{x}^H \mathbf{M}^{-1} \mathbf{x}} \right)^{-1}.$$

By defining $l_G(\mathbf{x}) = \frac{L_G(\mathbf{x})^{1/mL} - 1}{L_G(\mathbf{x})^{1/mL}}$ or $L_G(\mathbf{x})^{1/mL} = [1 - l_G(\mathbf{x})]^{-1}$, we obtain an equivalent test:

$$l_G(\mathbf{x}) = \frac{\mathbf{x}^H \mathbf{M}^{-1} \mathbf{P} (\mathbf{P}^H \mathbf{M}^{-1} \mathbf{P})^{-1} \mathbf{P}^H \mathbf{M}^{-1} \mathbf{x}}{\mathbf{x}^H \mathbf{M}^{-1} \mathbf{x}}, \quad (33)$$

which corresponds to the subspace version of the ACE test presented in [54].

## APPENDIX C

### RAO'S DETECTOR DERIVATION

The partial derivative of the log-likelihood function is defined as:

$$\frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \boldsymbol{\xi}_R, \boldsymbol{\xi}_S)}{\partial \boldsymbol{\xi}_R} = \left[ \frac{\frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \boldsymbol{\xi}_R, \boldsymbol{\xi}_S)}{\partial \text{Re}(\boldsymbol{\alpha})}}{\frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \boldsymbol{\xi}_R, \boldsymbol{\xi}_S)}{\partial \text{Im}(\boldsymbol{\alpha})}} \right].$$

From [53] (15.60), we obtain:

$$\frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \boldsymbol{\xi}_R, \boldsymbol{\xi}_S)}{\partial \boldsymbol{\xi}_R} = \left[ \begin{matrix} 2 \operatorname{Re}[\mathbf{P}^H \mathbf{C}^{-1}(\boldsymbol{\xi}_S) (\mathbf{x} - \mathbf{P}\boldsymbol{\alpha})] \\ 2 \operatorname{Im}[\mathbf{P}^H \mathbf{C}^{-1}(\boldsymbol{\xi}_S) (\mathbf{x} - \mathbf{P}\boldsymbol{\alpha})] \end{matrix} \right].$$

Thus:

$$\frac{\partial \ln p_{\mathbf{x}}(\mathbf{x}; \boldsymbol{\xi}_R, \boldsymbol{\xi}_S)}{\partial \boldsymbol{\xi}_R} \Big|_{\boldsymbol{\xi}_R = \mathbf{0}} = \left[ \begin{matrix} 2 \operatorname{Re}[\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{x}] \\ 2 \operatorname{Im}[\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{x}] \end{matrix} \right]. \quad (34)$$

Using [53] (15.52) $\mathbf{I}_{\boldsymbol{\xi}_S \boldsymbol{\xi}_R}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) = \mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) = \mathbf{0}$, thus $[\mathbf{I}^{-1}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S)]_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S} = \mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S)$. The elements of $\mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S)$ are given by:

$$[\mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S)]_{11} = 2 \mathbf{P}^H \mathbf{C}^{-1}(\boldsymbol{\xi}_S) \mathbf{P},$$

$$[\mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S)]_{22} = 2 \mathbf{P}^H \mathbf{C}^{-1}(\boldsymbol{\xi}_S) \mathbf{P},$$

$$[\mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S)]_{12} = 2 \operatorname{Re}[\mathbf{I} \mathbf{P}^H \mathbf{C}^{-1}(\boldsymbol{\xi}_S) \mathbf{P}] = \mathbf{0},$$

$$[\mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S)]_{21} = 2 \operatorname{Re}[-\mathbf{I} \mathbf{P}^H \mathbf{C}^{-1}(\boldsymbol{\xi}_S) \mathbf{P}] = \mathbf{0}.$$

We have therefore:

$$\mathbf{I}_{\boldsymbol{\xi}_R \boldsymbol{\xi}_S}(\boldsymbol{\xi}_R, \boldsymbol{\xi}_S) = \begin{bmatrix} 2 \mathbf{P}^H \mathbf{C}^{-1}(\boldsymbol{\xi}_S) \mathbf{P} & \mathbf{0} \\ \mathbf{0} & 2 \mathbf{P}^H \mathbf{C}^{-1}(\boldsymbol{\xi}_S) \mathbf{P} \end{bmatrix}. \quad (35)$$
Finally, by replacing (34) and (35) into (12) leads to:

$$L_R(\mathbf{x}) = \begin{bmatrix} 2 \operatorname{Re}(\mathbf{x}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{P}) & -2 \operatorname{Im}(\mathbf{x}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{P}) \\ \frac{(\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{P})^{-1}}{2} & \frac{\mathbf{0}}{(\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{P})^{-1}} \\ \frac{\mathbf{0}}{2} & \frac{2 \operatorname{Re}(\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{x})}{2 \operatorname{Im}(\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{x})} \end{bmatrix},$$

which simplifies to:

$$L_R(\mathbf{x}) = 2 \begin{bmatrix} \operatorname{Re}(\mathbf{x}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{P}) (\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{P})^{-1} \\ \operatorname{Re}(\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{x}) - \operatorname{Im}(\mathbf{x}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{P}) \\ (\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{P})^{-1} \operatorname{Im}(\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{x}) \end{bmatrix}.$$

Knowing that $\mathbf{P}^H \mathbf{C}^{-1}(\hat{\boldsymbol{\xi}}_{S_0}) \mathbf{P}$ is real and positive definite, it can be factorized and incorporated into the real and imaginary parts. After some algebraic manipulation, we obtain (14).

# APPENDIX D

## MAXIMUM LIKELIHOOD ESTIMATOR OF THE COVARIANCE MATRIX IN COMPOUND-GAUSSIAN CLUTTER

The likelihood of $\{\mathbf{x}_k\}_{k \in [1,K]}$ under $H_0$ can be rewritten as:

$$\begin{aligned} p_{\mathbf{x}}(\{\mathbf{x}_k\}_{k \in [1,K]}; \mathbf{M}, \mathbf{T}_k, H_0) \\ = \prod_{k=1}^{K} \frac{1}{\pi^{2m} |\tilde{\mathbf{C}}|} \exp(-\mathbf{x}_k^H \tilde{\mathbf{C}}^{-1} \mathbf{x}_k), \\ = \frac{1}{\pi^{2m K} |\mathbf{M}|^K} \prod_{k=1}^{K} \frac{1}{|\mathbf{T}_k|^2} \exp(-\mathbf{x}_k^H \mathbf{T}_k^{-1} \mathbf{M}^{-1} \mathbf{T}_k^{-1} \mathbf{x}_k), \end{aligned}$$

$$\text{where } \tilde{\mathbf{C}} = \mathbf{T}_k \mathbf{M} \mathbf{T}_k, \text{ and } \mathbf{T}_k = \begin{bmatrix} \sqrt{\tau_{1k}} & 0 \\ 0 & \sqrt{\tau_{2k}} \end{bmatrix} \otimes \mathbf{I}_m.$$

The log-likelihood can be written as:

$$\ln p_{\mathbf{x}}(\{\mathbf{x}_k\}_{k \in [1,K]}; \mathbf{M}, \mathbf{T}_k, H_0) = -2mK \ln \pi - K \ln |\mathbf{M}| + 2 \sum_{k=1}^{K} \ln |\mathbf{T}_k^{-1}| - \sum_{k=1}^{K} \mathbf{x}_k^H \mathbf{T}_k^{-1} \mathbf{M}^{-1} \mathbf{T}_k^{-1} \mathbf{x}_k. \quad (36)$$

According to [52] (82), the derivative with respect to $\mathbf{T}_k^{-1}$ is:

$$\frac{\partial \ln p_{\mathbf{x}}(\{\mathbf{x}_k\}_{k \in [1,K]}; \mathbf{M}, \mathbf{T}_k, H_0)}{\partial \mathbf{T}_k^{-1}} = 2 \mathbf{T}_k - 2 \operatorname{Re}(\mathbf{M}^{-1} \mathbf{T}_k^{-1} \mathbf{x}_k \mathbf{x}_k^H). \quad (37)$$

Following the same approach as in Appendix A, we obtain the minimum for $\mathbf{T}_k$ for a fixed $\mathbf{M}$:

$$\hat{\mathbf{T}}_k = \operatorname{Re} (\mathbf{M}^{-1} \hat{\mathbf{T}}_k^{-1} \mathbf{x}_k \mathbf{x}_k^H), \quad (38)$$

where

$$\hat{\mathbf{T}}_k = \begin{bmatrix} \sqrt{\hat{\tau}_{1k}} & 0 \\ 0 & \sqrt{\hat{\tau}_{2k}} \end{bmatrix} \otimes \mathbf{I}_m, \quad (39)$$

and

$$\hat{\tau}_{1k} = t_1 + \sqrt{\frac{t_1}{t_2}} t_{12}, \quad (40)$$

$$\hat{\tau}_{2k} = t_2 + \sqrt{\frac{t_2}{t_1}} t_{12}, \quad (41)$$

with

$$t_1 = \frac{1}{m} \mathbf{x}_{1,k}^H \mathbf{M}_{11}^{-1} \mathbf{x}_{1,k}, \quad (42)$$

$$t_2 = \frac{1}{m} \mathbf{x}_{2,k}^H \mathbf{M}_{22}^{-1} \mathbf{x}_{2,k}, \quad (43)$$

$$t_{12} = \frac{1}{m} \operatorname{Re} (\mathbf{x}_{1,k}^H \mathbf{M}_{12}^{-1} \mathbf{x}_{2,k}). \quad (44)$$

Replacing $\mathbf{T}_k$ by $\hat{\mathbf{T}}_k$ in (36) and deriving with respect to $\mathbf{M}^{-1}$ lead to:

$$\frac{\partial \ln p_{\mathbf{x}}(\{\mathbf{x}_k\}_{k \in [1,K]}; \mathbf{M}, \hat{\mathbf{T}}_k, H_0)}{\partial \mathbf{M}^{-1}} = K \mathbf{M} - \sum_{k=1}^{K} (\hat{\mathbf{T}}_k^{-1} \mathbf{x}_k) (\hat{\mathbf{T}}_k^{-1} \mathbf{x}_k)^H, \quad (45)$$

and the minimum in $\mathbf{M}$ is given by:

$$\hat{\mathbf{M}} = \frac{1}{K} \sum_{k=1}^{K} (\hat{\mathbf{T}}_k^{-1} \mathbf{x}_k) (\hat{\mathbf{T}}_k^{-1} \mathbf{x}_k)^H, \quad (46)$$

$$= \frac{1}{K} \sum_{k=1}^{K} \hat{\mathbf{T}}_k^{-1} \mathbf{x}_k \mathbf{x}_k^H \hat{\mathbf{T}}_k^{-1}. \quad (47)$$

The estimator (47) is independent of the textures. This could be shown by substituting $\mathbf{x}_k = \begin{bmatrix} \mathbf{x}_{1,k} \\ \mathbf{x}_{2,k} \end{bmatrix}$ by $\begin{bmatrix} \sqrt{\tau_{1k}} \mathbf{c}_{1,k} \\ \sqrt{\tau_{2k}} \mathbf{c}_{2,k} \end{bmatrix}$ in (39) and (47).

# REFERENCES

[1] I. Karoui, I. Quidu, and M. Legris, “Automatic sea-surface obstacle detection and tracking in forward-looking sonar image sequences,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 53, no. 8, pp. 4661–4669, 2015.

[2] I. Quidu, L. Jaulin, A. Bertholom, and Y. Dupas, “Robust multitarget tracking in forward-looking sonar image sequences using navigational data,” *IEEE Journal of Oceanic Engineering*, vol. 37, no. 3, pp. 417–430, 2012.

[3] I. Quidu, A. Bertholom, and Y. Dupas, “Ground obstacle tracking on forward looking sonar images,” in *Proceedings of ECUA’10*, July 2010.

[4] Y. Petillot, I. Tena Ruiz, and D. M. Lane, “Underwater vehicle obstacle avoidance and path planning using a multi-beam forward-looking sonar,” *IEEE Journal of Oceanic Engineering*, vol. 26, no. 2, pp. 240–251, 2001.
[5] B. Braginsky and H. Guterman, “Obstacle avoidance approaches for autonomous underwater vehicle: Simulation and experimental results,” *IEEE Journal of Oceanic Engineering*, vol. 41, no. 4, pp. 882–892, 2016.

[6] D. A. Abraham and P. K. Willett, “Active sonar detection in shallow water using the page test,” *IEEE Journal of Oceanic Engineering*, vol. 27, no. 1, pp. 35–46, 2002.

[7] D. A. Abraham, “Detection-threshold approximation for non-Gaussian backgrounds,” *IEEE Journal of Oceanic Engineering*, vol. 35, no. 2, pp. 355–365, 2010.

[8] D. A. Abraham, “Distributed active sonar detection in dependent *K*-distributed clutter,” *IEEE Journal of Oceanic Engineering*, vol. 34, no. 3, pp. 343–357, 2009.

[9] E. J. Kelly, “An adaptive detection algorithm,” *IEEE Transactions on Aerospace and Electronic Systems*, vol. AES-22, no. 2, pp. 115–127, 1986.

[10] F. C. Robey, D. R. Fuhrmann, E. J. Kelly, and R. Nitzberg, “A CFAR adaptive matched filter detector,” *IEEE Transactions on Aerospace and Electronic Systems*, vol. 28, no. 1, pp. 208–216, 1992.

[11] L. L. Scharf and B. Friedlander, “Matched subspace detectors,” *IEEE Transactions on Signal Processing*, vol. 42, no. 8, pp. 2146–2157, 1994.

[12] S. Kraut, L. L. Scharf, and L. T. McWhorter, “Adaptive subspace detectors,” *IEEE Transactions on Signal Processing*, vol. 49, no. 1, pp. 1–16, 2001.

[13] J. Ward, “Space-Time Adaptive Processing for airborne radar,” Tech. Rep., Lincoln Lab., MIT, Lexington, Mass., USA, December 1994.

[14] E. Fishler, A. Haimovich, R. S. Blum, L. J. Cimini, D. Chizhik, and R. A. Valenzuela, “Spatial diversity in radars—models and detection performance,” *IEEE Transactions on Signal Processing*, vol. 54, no. 3, pp. 823–838, 2006.

[15] E. Fishler, A. Haimovich, R. Blum, R. Cimini, D. Chizhik, and R. Valenzuela, “Performance of MIMO radar systems: advantages of angular diversity,” in *Conference Record of the Thirty-Eighth Asilomar Conference on Signals, Systems and Computers, 2004.*, 2004, vol. 1, pp. 305–309 Vol.1.

[16] N. H. Lehmann, A. M. Haimovich, R. S. Blum, and L. Cimini, “High resolution capabilities of MIMO radar,” in *2006 Fortieth Asilomar Conference on Signals, Systems and Computers*, 2006, pp. 25–30.

[17] J. Li and P. Stoica, “MIMO radar with colocated antennas,” *IEEE Signal Processing Magazine*, vol. 24, no. 5, pp. 106–114, 2007.

[18] L. Wei, W. Zhijie, L. Jianchen, W. Mingzhou, H. Qiao, and Y. Baomin, “Application of improved space-time adaptive processing in sonar,” *2012 International Conference on Computer Science and Electronics Engineering*, vol. 3, pp. 344–348, 2012.

[19] W. Li, G. Chen, E. Blasch, and R. Lynch, “Cognitive MIMO sonar based robust target detection for harbor and maritime surveillance applications,” in *2009 IEEE Aerospace conference*, 2009, pp. 1–9.

[20] N. M. Sasi, P. S. Sathidevi, R. Pradeepa, and S. Gopi, “A low complexity STAP for reverberation cancellation in active sonar detection,” in *2010 IEEE Sensor Array and Multichannel Signal Processing Workshop*, 2010, pp. 245–248.

[21] S. U. Pillai, J. R. Guerci, and S. R. Pillai, “Wideband STAP (WB-STAP) for passive sonar,” in *Oceans 2003. Celebrating the Past ... Teaming Toward the Future (IEEE Cat. No.03CH37492)*, 2003, vol. 5, pp. 2814–2818.

[22] X. Pan, Z. Liu, P. Zhang, Y. Shen, and J. Qiu, “Distributed MIMO sonar for detection of moving targets in shallow sea environments,” *Applied Acoustics*, vol. 185, pp. 108366, 2022.

[23] X. Liu, T. Wei, C. Sun, Y. Yang, and J. Zhuo, “High-resolution two-dimensional imaging using MIMO sonar with limited physical size,” *Applied Acoustics*, vol. 182, pp. 108280, 2021.

[24] J. D. Tucker and M. R. Azimi-Sadjadi, “Coherence-based underwater target detection from multiple disparate sonar platforms,” *IEEE Journal of Oceanic Engineering*, vol. 36, no. 1, pp. 37–51, 2011.

[25] J. Xue, M. Ma, J. Liu, M. Pan, S. Xu, and J. Fang, “Wald- and Rao-based detection for maritime radar targets in sea clutter with lognormal texture,” *IEEE Transactions on Geoscience and Remote Sensing*, vol. 60, pp. 1–9, 2022.

[26] D. Ciunzo, V. Carotenuto, and A. De Maio, “On multiple covariance equality testing with application to SAR change

detection,” *IEEE Transactions on Signal Processing*, vol. 65, no. 19, pp. 5078–5091, 2017.

[27] M. Sun, W. Liu, J. Liu, and C. Hao, “Complex parameter Rao, Wald, Gradient, and Durbin tests for multichannel signal detection,” *IEEE Transactions on Signal Processing*, vol. 70, pp. 117–131, 2022.

[28] P. Wang, H. Li, and B. Himed, “A parametric moving target detector for distributed MIMO radar in non-homogeneous environment,” *IEEE Transactions on Signal Processing*, vol. 61, no. 9, pp. 2282–2294, 2013.

[29] W. Liu, W. Xie, J. Liu, and Y. Wang, “Adaptive double subspace signal detection in Gaussian background—part I: Homogeneous environments,” *IEEE Transactions on Signal Processing*, vol. 62, no. 9, pp. 2345–2357, 2014.

[30] B. Y. Mills and A. G. Little, “A high-resolution aerial system of a new type,” *Australian Journal of Physics*, vol. 6, no. 3, pp. 272–278, 1953.

[31] E. Ollila, D. E. Tyler, V. Koivunen, and H. V. Poor, “Complex elliptically symmetric distributions: Survey, new results and applications,” *IEEE Transactions on Signal Processing*, vol. 60, no. 11, pp. 5597–5625, 2012.

[32] A. Farina, F. Gini, M.V. Greco, and L. Verrazzani, “High resolution sea clutter data: statistical analysis of recorded live data,” *IEE Proceedings - Radar, Sonar and Navigation*, vol. 144, pp. 121–130(9), June 1997.

[33] J. B. Billingsley, A. Farina, F. Gini, M. V. Greco, and L. Verrazzani, “Statistical analyses of measured radar ground clutter data,” *IEEE Transactions on Aerospace and Electronic Systems*, vol. 35, no. 2, pp. 579–593, 1999.

[34] M. Greco, F. Bordoni, and F. Gini, “X-band sea-clutter non-stationarity: influence of long waves,” *IEEE Journal of Oceanic Engineering*, vol. 29, no. 2, pp. 269–283, 2004.

[35] J.-P. Ovarlez, F. Pascal, and P. Forster, “Covariance Matrix Estimation in SIRV and Elliptical Processes and Their Applications in Radar Detection,” in *Modern Radar Detection Theory*, IET A. De Maio and M. S. Greco, Eds., number chapter 8, pp. 295–332. Scitech Publishing, 2015.

[36] A.-A. Saucan, T. Chonavel, C. Sintes, and J.-M. Le Caillec, “CPHD-DOA tracking of multiple extended sonar targets in impulsive environments,” *IEEE Transactions on Signal Processing*, vol. 64, no. 5, pp. 1147–1160, 2016.

[37] A. Hippert-Ferrer, M.N. El Korso, A. Breloy, and G. Ginolhac, “Robust low-rank covariance matrix estimation with a general pattern of missing values,” *Signal Processing*, vol. 195, pp. 108460, 2022.

[38] F. Pascal, Y. Chitour, J.-P. Ovarlez, P. Forster, and P. Larzabal, “Covariance structure maximum-likelihood estimates in compound Gaussian noise: Existence and algorithm analysis,” *IEEE Transactions on Signal Processing*, vol. 56, no. 1, pp. 34–48, 2008.

[39] M. Mahot, F. Pascal, P. Forster, and J.-P. Ovarlez, “Asymp totic properties of robust complex covariance matrix estimates,” *IEEE Transactions on Signal Processing*, vol. 61, no. 13, pp. 3348–3356, 2013.

[40] C. Y. Chong, F. Pascal, J.-P. Ovarlez, and M. Lesturgie, “MIMO radar detection in non-Gaussian and heterogeneous clutter,” *IEEE Journal of Selected Topics in Signal Processing*, vol. 4, pp. 115–126, March 2010.

[41] N. Petrov, F. Le Chevalier, and A. G. Yarovoy, “Detection of range migrating targets in compound-Gaussian clutter,” *IEEE Transactions on Aerospace and Electronic Systems*, vol. 54, no. 1, pp. 37–50, 2018.

[42] F. Mosca, G. Matte, O. Lerda, F. Naud, D. Charlot, M. Rioblanc, and C. Corbières, “Scientific potential of a new 3D multibeam echosounder in fisheries and ecosystem research,” *Fisheries Research*, vol. 178, pp. 130–141, 2016.

[43] C. Corbières and F. Mosca, “Real-time acoustic fishing selectivity for sustainable practices,” in *OCEANS 2017 - Anchorage*, 2017, pp. 1–8.

[44] B. Tallon, P. Roux, G. Matte, J. Guillard, J. H. Page, and E. Skipetrov, “Ultra slow acoustic energy transport in dense fish aggregates,” *Scientific Reports*, vol. 11, no. 17541, September 2021.

[45] G. Matte, D. Charlot, O. Lerda, T.-K. N’Guyen, V. Giovanini, M. Rioblanc, and F. Mosca, “Seapix: an innovative multibeam multiswath echosounder for water column and seabed analysis,” in *Hydro17*, 2017.
[46] T.-K. Nguyen, D. Charlot, J.-M. Boucher, G. Le Chenadec, and R. Fablet, "Seabed classification using a steerable multibeam echo sounder," in *OCEANS 2016 MTS/IEEE Monterey*, 2016, pp. 1-4.

[47] A. Aubry, J. Carretero-Moya, A. De Maio, A. Pauciullo, J. Gismero-Menoyo, and A. Asensio-Lopez, "Detection of Extended Target in Compound-Gaussian Clutter," in *Modern Radar Detection Theory*, IET A. De Maio and M. S. Greco, Eds., number chapter 9, pp. 333-374. Scitech Publishing, 2016.

[48] S. M. Kay, *Fundamentals of Statistical Signal Processing: Detection theory*, Prentice-Hall PTR, 1998.

[49] S. Kay and Z. Zhu, "The complex parameter Rao test," *IEEE Transactions on Signal Processing*, vol. 64, no. 24, pp. 6580-6588, 2016.

[50] C. Y. Chong, F. Pascal, J.-P. Ovarlez, and M. Lesturgie, "Robust MIMO radar detection for correlated subarrays," in *IEEE International Conference on Acoustic, Speech and Signal Processing (ICASSP)*, April 2010, pp. 2786-2789.

[51] I. S. Reed, J. D. Mallett, and L. E. Brennan, "Rapid convergence rate in adaptive arrays," *IEEE Transactions on Aerospace and Electronic Systems*, vol. AES-10, no. 6, pp. 853-863, 1974.

[52] K. B. Petersen and M. S. Pedersen, *The Matrix Cookbook*, Technical University of Denmark, nov 2012.

[53] S. M. Kay, *Fundamentals of Statistical Signal Processing: Estimation Theory*, Prentice-Hall PTR, 1993.

[54] R. Raghavan, S. Kraut, and C. Richmond, *Subspace Detection for Adaptive Radar: Detectors and Performance Analysis*, pp. 43-83, Institution of Engineering and Technology, 2016.