IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO.

21, NOVEMBER 1, 2017 5649

Robust Calibration of Radio Interferometers in

Non-Gaussian Environment
Virginie Ollier, Mohammed Nabil El Korso, Rémy Boyer, Senior Member, IEEE, Pascal Larzabal, Member, IEEE,
and Marius Pesavento, Member, IEEE

Abstract—The development of new phased-array systems in ra- Index Terms—Calibration, robustness, spherically invariant
dio astronomy, as the low-frequency array (LOFAR) and the square random process, relaxed concentrated maximum likelihood, Jones
kilometre array (SKA), formed of a large number of small and flex- matrices, radio astronomy.
ible elementary antennas has led to significant challenges. Among
them, model calibration is a crucial step in order to provide ac-
curate and, thus, meaningful images and requires the estimation I. INTRODUCTION
of all the perturbation effects introduced along the signal propa- ADIO astronomy aims to study radio emissions from the
gation path for a specific source direction and antenna position.
Usually, it is common to perform model calibration using the a
priori knowledge regarding a small number of known strong cal-
R sky, in order to detect, identify new objects and observe
known structures at higher resolution, in a specific electromag-
ibrator sources but under the assumption of Gaussianity of the netic spectrum [1]. This fundamental thematic shines a new light
noise. Nevertheless, observations in the context of radio astronomy on our universe, revealing more about its nature and history. In
are known to be affected by the presence of outliers, which are order to carry out particularly sensitive observations in a large
due to several causes, e.g., weak non-calibrator sources or man-
made radio frequency interferences. Consequently, the classical range of the spectrum, and to handle significant cosmological
Gaussian noise assumption is violated leading to severe degrada- issues, largely distributed sensor arrays are currently being built
tion in performances. In order to take into account the outlier or planned, such as the low frequency array (LOFAR) [2] and
effects, we assume that the noise follows a spherically invariant the square kilometre array (SKA) [3]. They will notably be com-
random distribution. Based on this modeling, a robust calibra- posed of a large number of relatively low-cost small antennas
tion algorithm is presented in this paper. More precisely, this new
scheme is based on the design of an iterative relaxed concentrated with wide field of view, resulting in a large collecting area and
maximum likelihood estimation procedure that allows us to obtain high resolution imaging. Nevertheless, to meet the theoretical
closed-form expressions for the unknown parameters with a rea- optimal performances of such next generation radio interfer-
sonable computational cost. Numerical simulations reveal that the ometers, a plethora of signal processing challenges must be
proposed algorithm outperforms the state-of-the-art calibration overcome, among them, calibration, data reduction and image
techniques.
synthesis [4]–[7]. These aspects are intertwined and must be
dealt with to take advantage of the new advanced radio interfer-
ometers. As an example, lack of calibration has dramatic effects
Manuscript received December 6, 2016; revised May 21, 2017; accepted July in the image reconstruction by causing severe distortions. In this
10, 2017. Date of publication July 31, 2017; date of current version August
31, 2017. The associate editor coordinating the review of this manuscript and paper, we focus on calibration, which involves the estimation of
approving it for publication was Prof. Adel Belouchrani. This work was sup- all unknown perturbation effects and represents a cornerstone
ported in part by MAGELLAN (ANR-14-CE23-0004-01), in part by ON FIRE of the imaging step [8]–[10].
Project (Jeunes Chercheurs GDR-ISIS), and in part by ANR ASTRID Project
MARGARITA. (Corresponding author: Virginie Ollier.) Array calibration aspects have been tackled for a few decades
V. Ollier is with the Laboratory Systèmes et Applications des Technolo- in the array processing community leading to a variety of cal-
gies de l’Information et de l’Energie, UMR 8029, École Normale Supérieure ibration algorithms [11]–[13]. Such algorithms can be classi-
de Cachan, Université Paris-Saclay, Cachan 94235, France, and also with the
Laboratoire des Signaux et Systèmes UMR 8506, Université Paris-Sud, Univer- fied into two different approaches depending on the presence
sité Paris-Saclay, Gif-sur-Yvette 91192, France (e-mail: virginie.ollier@satie. [14]–[16], or the absence [17]–[22], of one or more cooperative
ens-cachan.fr). sources, named calibrator sources. In the radio astronomy con-
M. N. El Korso is with the LEME, EA 4416, Université Paris-Ouest, Ville
d’Avray 92410, France (e-mail: m.elkorso@parisnanterre.fr). text, calibration is commonly treated using the first approach as
R. Boyer is with the Laboratoire des Signaux et Systèmes, UMR 8506, we have access to prior knowledge thanks to tables describing
Université Paris-Sud, Université Paris-Saclay, Gif-sur-Yvette 91192, France (e- accuratly the position and flux of the brightest sources [23].
mail: remy.boyer@l2s.centralesupelec.fr).
P. Larzabal is with the Laboratory Systèmes et Applications des Technolo- Following this methodology, the majority of proposed cali-
gies de l’Information et de l’Energie, UMR 8029, École Normale Supérieure bration schemes in radio interferometry are least squares-based
de Cachan, Université Paris-Saclay, Cachan 94235, France (e-mail: pascal. approaches. The state-of-the-art consists in the so-called alter-
larzabal@satie.ens-cachan.fr).
M. Pesavento is with the Communication Systems Group, Technische nating least squares approach [24]–[27], which leads to statis-
Universität Darmstadt, Darmstadt 64289, Germany (e-mail: mpesa@nt.tu- tically efficient algorithm under a Gaussian model, since the
darmstadt.de). least squares estimator is equivalent to the maximum likelihood
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org. (ML) estimator in this case. On the other hand, expectation
Digital Object Identifier 10.1109/TSP.2017.2733496 maximization (EM) [28]–[30] and EM-based algorithms, such
1053-587X © 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications standards/publications/rights/index.html for more information.
5650 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 21, NOVEMBER 1, 2017

as the space alternating generalized expectation maximization effect separately thanks to individual Jones terms in a cumulative
algorithm [31], have been proposed in order to enhance the con- product. Thus, different corruptions are described by different
vergence rate of the least squares-based calibration algorithms kinds of Jones matrices. We emphasize that the proposed algo-
[32]. Nevertheless, the major drawback of these schemes is the rithm, entitled relaxed concentrated ML estimator, is a generic
Gaussianity assumption which is not realistic in the radio astron- algorithm as it is based on a non-structured Jones matrices for-
omy context. Specifically, the presence of outliers has multiple mulation as a first step. However, it can be adapted to various
causes, among which i) the radio frequency interferers, which regimes describing distinct calibration scenarios in which an
corrupt the observations and are not always perfectly filtered in array can operate [46]. In this paper, we consider the specific
practice [33], [34], ii) the presence of unknown weak sources in example of the direction dependent distortion regime with a
the background [35], iii) the presence of some punctual events compact set of antennas, which we refer to as the 3DC regime.
such as interference due to the Sun or due to strong sources in The array is therefore considered as a closely packed group of
the sidelobes which can also randomly create outliers [36]. To antennas but the array elements have a wide field of view. This
the best of our knowledge, the proposed scheme in [35], repre- is particularly well-adapted for calibration of compact arrays,
sents the only alternative to the existing calibration algorithms typically a LOFAR station.
based on a Gaussian noise model. The rest of the paper is organized as follows: in Section II,
In [35], theoretical and experimental analyses have been con- we present the data model in the context of radio astronomy,
ducted in order to demonstrate that the effect of outliers in first with non-structured Jones matrices and thereafter, we study
the radio astronomy context can indeed be modeled by a non- an example of structured Jones matrices for the 3DC calibra-
Gaussian heavy-tailed distributed noise process. Nevertheless, tion regime. In Section III, we give an overview of the proposed
the algorithm presented in [35] has its own limits, since the robust ML estimator, based on spherically invariant random pro-
noise is specifically modeled as a Student’s t with independent cess (SIRP) noise modeling. An efficient estimation procedure
identically distributed entries. To improve the robustness of the of the distortions introduced on each signal propagation path
calibration, we propose, in this paper, a new scheme based on is derived in Section IV. Then, the algorithm is adapted to the
a broader class of distributions gathered under the so-called case of structured Jones matrices in Section V for the 3DC cal-
spherically invariant random noise modeling [37], [38], which ibration regime. Finally, we provide numerical simulations in
includes the Student’s t distribution. A spherically invariant ran- Section VI to assess the robustness of the approach and draw
dom vector (SIRV) is described as the product of a positive our conclusions in Section VII.
random variable, named texture, and the so-called speckle com- In this paper, we use the following notation: symbols (·)T ,
ponent which is Gaussian, resulting in a two-scale compound (·)∗ , (·)H denote, respectively, the transpose, the complex conju-
Gaussian distribution [39]. The flexibility of the SIRV modeling gate and the Hermitian transpose. The Kronecker product is rep-
allows to consider non-Gaussian heavy-tailed distributed noise resented by ⊗, E{·} denotes the expectation operator, bdiag{·}
in the presence of outliers, but also to adaptively consider Gaus- is the block-diagonal operator, whereas diag{·} converts a vec-
sian noise in the extreme case when there are no outliers. Under tor into a diagonal matrix. The trace and determinant operators
the SIRV model, we estimate the unknown parameters itera- are, respectively, referred by tr {·} and | · |. The symbol IB rep-
tively based on a relaxed ML estimator, leading to closed-form resents the B × B identity matrix, vec(·) stacks the columns of
expressions for the noise parameters while a block coordinate a matrix on top of one another, || · ||F is the Frobenius norm,
descent (BCD) algorithm [40], [41] is designed to obtain the while || · ||2 denotes the l2 norm. Finally,  {·} represents the
estimates of parameters of interest efficiently and at a low cost. real part and we note j the complex number whose square
Finally, it is worth mentioning that the parametric model used equals −1.
in this paper to describe the perturbation effects is based on the
so-called Jones matrices [42], [43]. Such formalism describes II. DATA MODEL
in a flexible way the conversion of the incident electric field into
voltages. Indeed, along its propagation path, the signal is af- A. Case of Non-Structured Jones Matrices
fected by various effects and transformations which correspond Let us consider M antennas with known locations that receive
to matrix multiplications in the mathematical Jones framework. D signals emitted by calibrator sources. Each antenna is dual
Multiple distortion effects caused by the environment and/or the polarized and composed of two receptors, in order to provide
instruments can be easily incorporated into the model using an sensitivity to the two orthogonal polarization directions (x, y)
adequate parametrization of the Jones matrices. Such effects can of the incident electromagnetic plane wave. Consequently, the
represent, for example, the ionospheric phase delay resulting in relation between the i-th source emission and the measured
angular shifts, the atmospheric distortions, the typical phase de- voltage at the p-th antenna is given by [42], [43], [47]
lay due to geometric pathlength difference, the voltage primary
vi,p (θ) = Ji,p (θ)si (1)
beam, the cross-leakage or also the electronic gains [44], [45].
T
For the above reasons and due to its flexibility [1], [32], [42]– where si = [si x , si y ] is the incoming signal, vi,p (θ) = [vi,p x
[44], we adopt this parametric model. We make a distinction (θ), vi,p y (θ)]T is the generated voltage with one output for
between the non-structured and the structured cases: in the first each polarization direction and Ji,p (θ) denotes the so-called
one, one total Jones matrix stands for all the effects along the full 2 × 2 Jones matrix, parametrized by the unknown vector of
signal path while in the second case, we regard each physical interest θ. The Jones matrix models the array response and all
OLLIER et al.: ROBUST CALIBRATION OF RADIO INTERFEROMETERS IN NON-GAUSSIAN ENVIRONMENT 5651

the perturbations introduced along the path from the i-th source zero-mean Gaussian distribution,1 i.e.,
to the p-th sensor. Since each propagation path is particular, we
can associate a different Jones matrix with each source-antenna gpq ∼ CN (0, Ω). (6)
pair (i, p), leading to a total number of DM Jones matrices. In order to remove scaling ambiguities, we impose tr {Ω} = 1.
In this section, we consider the non-structured case where no Note that the choice of this constraint is arbitrary and does not
specific perturbation model is used to describe the physical affect the estimates of interest as argued in [49].
mechanism behind each perturbation effect and the unknown In this section, we adopted the non-structured Jones matrices
elements correspond to the entries of all Jones matrices [32], formulation which is relevant in the radio astronomical context
[48] (a structured example is given for 3DC calibration regime, [32], [48]. In this case, there is no need to specify the full prop-
in Section II-B). agation path, avoiding misspecification in the model. Besides, it
For each antenna pair, we compute the correlation of the is highly flexible and can be adapted to different scenarios [46].
output signals, resulting in the typical observations recorded In the following, we present the direction dependent distortion
by a radio interferometer. The correlation between voltages is regime with a compact set of antennas, named 3DC regime.
given, in the case of noise free measurements, for the (p, q)
antenna pair, by B. Specific Case of the 3DC Calibration Regime
 D
 D

  For a specific propagation path, from the i-th source to the
H
Vpq (θ) = E vi,p (θ) vi,q (θ)
p-th antenna, the global Jones matrix Ji,p accounts for mul-
i=1 i=1
tiple effects which can be described explicitly. Indeed, each
D
 global matrix can be decomposed into individual Jones terms
= Ji,p (θ)Ci JH
i,q (θ) for p < q, p, q ∈ {1, . . . , M }, which stand for specific physical effects [43], [44]. This way,
i=1
(2) instead of estimating entries of all Jones matrices as done in
the non-structured case, we will estimate physically meaning-
where the signals emitted by the sources are assumed uncorre- ful parameters, thus reducing the total number of parameters to
lated and the 2 × 2 matrix Ci = E{si sHi } is known from prior estimate. Introducing structured Jones matrices can be done in
knowledge. Let us remark that autocorrelations are not consid- the context of calibration scenarios [10], [46]. In what follows,
ered as shown by the condition p < q in (2) (this is a typical we target one specific commonly used calibration regime that
situation in the radio astronomy context where radio interfero- we call 3DC calibration regime, described in Fig. 1, which is
metric systems automatically flag the autocorrelations [2]). well adapted for calibration of LOFAR on station level [26], and
Using the property vec(ABC) = (CT ⊗ A)vec(B), we also for stations of the future SKA radio interferometer [10]. In
rewrite (2) as a 4 × 1 vector this scenario, direction dependent distortions play a significant
  D role since individual elements in the array have a wide field
ṽpq (θ) = vec Vpq (θ) = ui,pq (θ) (3) of view. Indeed, this implies different propagation conditions
i=1 towards distinct sources in the field of view. However, the ar-
in which ui,pq (θ) = (J∗i,q (θ) ⊗ Ji,p (θ))ci with ci = vec(Ci ). ray being relatively compact, made of similar elements, some
We stack all the noisy measurements within a full vector x = effects might be the same for all antennas. In the following, we
introduce a particular sequence of Jones matrices with specific
T
[v12 T
, v13 T
, . . . , v(M T
−1)M ] ∈ C
4B ×1
, where B = M (M2 −1) de-
parametrizations, in the context of 3DC calibration regime [50]
notes the total number of antenna pairs and vpq = ṽpq (θ) + npq
with npq the noise sample at a specific antenna pair. Specifically, i,p ) = Gp (gp )Hi,p Zi,p (αi )Fi (ϑi )
Ji,p (θ 3DC (7)
x reads
T
D
for i ∈ {1, . . . , D}, p ∈ {1, . . . , M } and θ 3DC
i,p = [ϑi , gp ,

x= ui (θ) + n (4) αTi ]T . We note Hi,p the only assumed known matrix thanks
i=1 to electromagnetic simulations in terms of antenna response
and a priori knowledge given by calibrator source and antenna
in which ui (θ) = [uTi,12 (θ), uTi,13 (θ), . . . , uTi,(M −1)M (θ)]T positions [43], [44], [50], [51], whereas the remaining matrices
and n = [nT12 , nT13 , . . . , nT(M −1)M ]T is the full noise vector are explained in the following items.
which accounts for Gaussian noise, but also the presence of • Ionospheric effects: Propagation through the ionosphere,
outliers in our data. Therefore, the noise can no longer be con- the outer layer of the earth’s atmosphere, introduces propaga-
sidered Gaussian and a robust calibration method is required. tion delay on the signal which is affected by spatially variable
To investigate non-Gaussian noise modeling and encompass a fluctuations. If these perturbation effects are not corrected for,
broad range of noise distributions, we propose to adopt the SIRP the sources may appear shifted from their intrinsic positions
noise model [37], [38]. Specifically, the noise at each antenna [10], [52]. In the case of a compact array, the ionospheric delay
pair is assumed to be generated as
√
npq = τpq gpq , (5) 1 Let us note that it is possible to consider a different covariance matrix Ω
pq
for each speckle component g p q in (6). In this case, the proposed algorithm
where the positive real random variable τpq is referred to as requires a few modifications and the corresponding expressions are presented
texture, whereas the complex speckle component gpq follows a in Appendix A.
5652 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 21, NOVEMBER 1, 2017

III. ROBUST CALIBRATION ESTIMATOR

This section is devoted to the design of a robust calibration
estimator based on the model (4). As it can be seen from (5),
one has to specify the probability density function (pdf) of each
texture parameter τpq in order to obtain the exact ML estimates.
Nevertheless, in pratical scenarios, such prior knowledge is not
available. Consequently, our idea is to make use of a relaxed
version of the exact model, i.e., we assume deterministic but un-
known texture realizations in the estimation process [53], [54].
This ensures more flexibility in our algorithm as the texture
distribution is not precisely described and avoids any possible
model misspecification, which is consistent with our motivation
to design a broad robust estimator w.r.t. the presence of outliers.
On the other hand, we adopt here an iterative procedure in or-
der to reduce the computational cost. In doing so, the proposed
algorithm sequentially updates each block of unknown param-
eters while fixing the remaining parameters. This leads to a
relaxed concentrated ML based calibration estimator for which
the expression of the likelihood function, when independency
is assumed between measurements, is written as
 
1 1 H −1
f (x|θ, τ , Ω) = exp − a (θ)Ω apq (θ) ,
pq
|πτpq Ω| τpq pq
(10)
Fig. 1. 3DC calibration regime, for which V  S and S  A. All receiving
elements in the station see the same ionosphere part but, due to their wide field where the vector composed of all texture realizations is τ =
of view, a multitude of sources are visible and perturbations are highly direction [τ12 , τ13 , . . . , τ(M −1)M ]T and apq (θ) = vpq − ṽpq (θ). Conse-
dependent. quently, the log-likelihood function reads

log f (x|θ, τ , Ω) = −4B log π

matrix is in fact a scalar direction-dependent phase given by
  1

−4 log τpq − B log |Ω| − aH (θ)Ω−1 apq (θ).
Zi,p (αi ) = exp jϕi,p I2 (8) pq pq
τpq pq
(11)
in which ϕi,p = ηi up + ζi vp where αi = [ηi , ζi ]T is the vector
of unknown offsets resulting in a shift of the i-th source direction In the following, we present the sequential updates of each
and rp = [up , vp ]T is the vector of known antenna position in block of unknown parameters, namely, τ , Ω and θ, following
units of wavelength. the methodology as in [55], [56].
On top of that, passing through the ionosphere is associated 1) Derivation of τ̂pq : Taking the derivative of the log-
with a rotation of the polarization plane of each signal source likelihood function in (11) w.r.t. τpq leads to the following tex-
around the line of sight. We call it the ionospheric Faraday ture estimate
rotation matrix Fi (ϑi ) and write it as 1 H
τ̂pq = a (θ)Ω−1 apq (θ). (12)
  4 pq
cos(ϑi ) − sin(ϑi )
Fi (ϑi ) = (9) 2) Derivation of Ω̂: The derivative of the log-likelihood func-
sin(ϑi ) cos(ϑi )
tion w.r.t. the element [Ω]k ,l of the speckle covariance matrix,
using classical differential properties [57, p. 2741], leads to
where ϑi is the unknown Faraday rotation angle, assumed iden-
tical for all antennas, since the array has a limited spatial extent  −1   1 −1 −1
−Btr Ω̂ ek eTl + aH T
pq (θ)Ω̂ ek el Ω̂ apq (θ) = 0
[44].
pq
τ pq
• Instrumental effects: Individual antennas are described by (13)
electronic complex gains which appear in Gp (gp ) = diag{gp } where the vector ek contains zeros except at the k-th position
with gp the unknown electronic gain vector. which is equal to unity. The permutation property of the trace
Therefore, in this specific structured case, the physi- operator enables to rewrite (13) as
cal model parameters in (7) are collected in the vector
T
3DC T 3DC T T −1  1 −1 −1
ε3DC = P[θ 3DC 1,1 , θ 1,2 , . . . , θ D ,M ] where P is an appro- −BeTl Ω̂ ek + eTl Ω̂ apq (θ)aH
pq (θ)Ω̂ ek = 0.
priate rearrangement matrix such that ε3DC = [ϑ1 , . . . , ϑD , g1T , pq
τ pq
T
. . . , gM , αT1 , . . . , αTD ]T . (14)
OLLIER et al.: ROBUST CALIBRATION OF RADIO INTERFEROMETERS IN NON-GAUSSIAN ENVIRONMENT 5653

This finally leads to the following estimate of the speckle co- function of this conditional distribution is maximized in the
variance matrix M-step. Therefore, this last step consists in an optimization pro-
1  1 cess which can be performed numerically, for instance with
Ω̂ = apq (θ)aH
pq (θ). (15) the Levenberg-Marquardt (LM) algorithm [59]–[61], or analyt-
B pq τpq
ically if closed-form expressions are available. As we show in
Inserting (12) into (15), we obtain the following, the optimization step is carried out w.r.t. to θ i
∈ C 4M ×1 instead of θ ∈ C 4D M ×1 . Therefore, the global mul-
h+1 4  apq (θ)aH
pq (θ)
Ω̂ =  h −1 (16) tiple source estimation problem is reduced to multiple single
B pq
aH
pq (θ) Ω̂ apq (θ) source sub-problems.
1) E-step: For the i-th source, we introduce the so-called
where h denotes the h-th iteration. Due to the introduced con- complete data vector wi such that
straint, we normalize the estimate of Ω by its trace, as follows
D

h+1
h+1 Ω̂ x= wi (20)
Ω̂ =  h+1  . (17) i=1
tr Ω̂

with wi = ui (θ i ) + ni and n = D i=1 ni , in which ni ∼
3) Estimation of θ̂: For given Ω and τ , estimating θ̂ is equiv- D
CN (0, βi Ψ). We have i=1 βi = 1 and Ψ is the covariance
alent to the following minimization problem matrix of n. Since npq ∼ CN (0, τpq Ω) and with the indepen-
 
 1 dence property, we obtain the following block-diagonal expres-
H −1
θ̂ = argmin a (θ)Ω apq (θ) . (18) sion for Ψ
θ pq
τpq pq
Ψ = bdiag{τ12 Ω, . . . , τ(M −1)M Ω}. (21)
In the following, we aim to reduce the computational cost of the
minimization procedure in (18) by use of the EM algorithm. For Let us note w = [w1T , . . . , wD
T T
] the complete data vector,
generality, we first adopt the non-structured Jones matrix formu- whose covariance matrix, denoted as Ξ, has the following form
lation, which can also be specified depending on the scenario,
as shown in Section V. Ξ = bdiag{β1 Ψ, . . . , βD Ψ}. (22)

IV. ESTIMATION OF θ̂ FOR NON-STRUCTURED With [62, p. 36], and after some calculus, the expression of
JONES MATRICES the conditional expectation is given by
Estimating directly the entries of the Jones matrices avoids ŵi = E{wi |x; θ, τ , Ω}
specifying any particular physical model, leading to a calibration  D

algorithm which is less sensitive to model errors in comparison 
= ui (θ i ) + βi x− ul (θ l ) . (23)
with algorithms based on the structured case. Due to the possible
l=1
large size of θ, a multi-dimensional parameter search needs to be
carried out to solve the optimization problem in (18) which re- 2) M-step: The goal of this step is to estimate θ i . Once ŵi
quires significant computation time. To reduce this complexity, are computed for i ∈ {1, . . . , D}, the estimated complete data
we make use of the EM algorithm. vector ŵ can be evaluated. The M-step is an optimization prob-
The essence of this algorithm relies on a proper parameter lem based on the following likelihood function where wi are
vector partitioning as well as an adequate choice of the so-called independent
complete data. As mentioned above, the parameters of interest
θ represent the entries of all Jones matrices. Consequently, it is f (ŵ |θ, τ , Ω)
natural to consider the following partition   H  
1
= exp − ŵ − u(θ) Ξ−1 ŵ − u(θ) (24)
θ = [θ T1 , . . . , θ TD ]T = [θ T1,1 , . . . , θ T1,M , . . . , θ TD ,1 , . . . , θ TD ,M ]T , |πΞ|
(19) D
  H  
1
for which the vector θ i,p ∈ R8×1 is the parametrization of the = exp − ŵi − ui (θ i ) (βi Ψ)−1 ŵi − ui (θ i ) .
path from the i-th calibrator source to the p-th sensor, i.e., i=1
|πβi Ψ|
Ji,p (θ) = Ji,p (θ i,p ).
To obtain an estimation of θ i , we need to minimize the following
cost function
A. Use of the EM Algorithm to Solve (18)
 H  
The EM algorithm [28]–[30] enables to compute the ML es- φi (θ i ) = ŵi − ui (θ i ) (βi Ψ)−1 ŵi − ui (θ i ) . (25)
timates and reduce the computational cost, via the iteration of
two steps. The first one is the E-step which reduces, in our To decrease even more the complexity cost of the proposed
scenario, to the computation of the conditional expectation of robust calibration scheme, we use the BCD algorithm [40], [41]
the complete data given the observed data and the current es- in the M-step. Consequently, we obtain analytical solutions for
timate of parameters [32], [58]. Afterwards, the log-likelihood each single source sub-problems in (25), as shown below.
5654 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 21, NOVEMBER 1, 2017

B. Use of the BCD Algorithm to Minimize (25)

In (25), the optimization is performed w.r.t. θ i . However, this
unknown parameter vector can be partitioned according to the
antennas, as expressed in (19). In the following, we perform the
optimization of the cost function w.r.t. each θ i,p (p-th antenna),
given the current estimates of all the other θ i,q with q = p.
This leads to a closed-form expression of θ̂ i,p as function of
θ i,q for q = p and the optimization process is repeated for each
component vector θ i,p for p ∈ {1, . . . , M }.
If we want to minimize (25) w.r.t. block-coordinate vector
θ i,p , we notice that only a subset of ui (θ i ) is actually depen-
dent on θ i,p , i.e., {ui,pq } for q > p, q ∈ {1, . . . , M } and
{ui,q p } for q < p, q ∈ {1, . . . , M }. Therefore, (25) reads
φi (θ i,p )
M 
 H  
= wi,pq − ui,pq (θ i,p ) (βi τpq Ω)−1 wi,pq − ui,pq (θ i,p )
q =1
q>p

M 
 H  
+ wi,q p − ui,q p (θ i,p ) (βi τq p Ω)−1 wi,q p − ui,q p (θ i,p ) iterations is considered in each loop, so we might not attain local
q =1
q<p convergence. However, we show in Section VI-A the relatively
good numerical stability of the algorithm.
+ Constant. (26) The scheme of the proposed algorithm is described in
By Constant, we mean the expressions independent of θ i,p . We Algorithm 1.
show in Appendix B that it is possible to write ui,pq directly as
a function of θ i,p , i.e.,
V. STRUCTURED JONES MATRICES
ui,pq (θ i,p ) = Σi,q θ i,p . (27)
We recall that the output of Algorithm 1 is the estimate of
In the same way, we have each Jones matrix denoted by Ĵi,p for i ∈ {1, . . . , D} and p ∈
ui,q p (θ i,p ) = Υi,q θ ∗i,p . (28) {1, . . . , M }. In the following, we consider the data model in (7)
for the specific 3DC calibration regime, and intend to estimate
Notation and calculations being introduced in Appendix B, we the unknown parameter vector of interest ε3DC in a sequential
only present here the results obtained, i.e., the expression of the manner. To do so, we use an iterative estimation procedure by
estimated entries of the Jones matrix associated with the path optimizing a cost function w.r.t. one parameter while fixing the
from the i-th calibrator source to the p-th sensor which is given others.
by 1) Estimation of gp : The diagonal elements of the gain matrix
⎧ are given by
⎪
⎪ (ΣH H
i Ai,p Σi + Υi Ãi,p Υi )
−1
×
⎪
⎪
⎨ (ΣH A w + ΥH Ã w̃ ) for 1 < p < M
i i,p i,p i i,p i,p
θ̂ i,p = ĝp = argmin κ(gp ) (30)
⎪
⎪ (Σ H
A Σ
i,p i )−1 H
Σ A i,p wi,p for p = 1 gp
⎪
⎪ i i
⎩ H H
(Υi Ãi,p Υi )−1 Υi Ãi,p w̃i,p for p = M
(29) 
where κ(gp) = D i=1 ||Ĵi,p −Gp(gp )Hi,p Zi,p Fi ||F . We rewrite
2
Therefore, θ i,p for p ∈ {1, . . . , M } are estimated in an iterative the cost function as
loop. With (25) and (29), it can be proven that the BCD algo-
rithm leads to unique solutions and thus, convergence to at least 
D
  H 
a local maximizer, is ensured [63]. If the M-step is performed κ(gp ) = Tr Ĵi,p − Gp (gp )Ri,p Ĵi,p − Gp (gp )Ri,p
exactly (i.e., the BCD gives the exact minimizer of (25) and con- i=1
sequently, the M step is exactly solved), convergence of the EM (31)
algorithm to a stationary point is ensured (to avoid the unusual in which Ri,p = Hi,p Zi,p Fi . The derivation of κ(gp ) w.r.t. the
case of convergence to a saddle point, a proper initialization is k-th element [gp ]k leads to
required) [30], with a theoretical infinite number of iterations.
Finally, in this case, convergence of the concentrated MLE is
∂κ(gp )   
D
guaranteed for an infinite number of iterations since the value = Tr − ek eTk Ri,p ĴH
i,p + ek eT
k Ri,p RH
G
i,p p
H
.
of the cost function at each step can either improve or main- ∂[gp ]k i=1
tain but cannot worsen [64]. In practice, only a finite number of (32)
OLLIER et al.: ROBUST CALIBRATION OF RADIO INTERFEROMETERS IN NON-GAUSSIAN ENVIRONMENT 5655

Let us denote Xi,p = Ri,p ĴH H

i,p and Wi,p = Ri,p Ri,p . From
(32), we deduce the equation satisfied by the gain matrix
D
 D

   
Xi,p k ,k
= Wi,p ĜH
p k ,k
i=1 i=1
D
    
= Wi,p k ,k
Ĝ∗p k ,k
(33)
i=1

for k ∈ {1, 2}, since Gp is a diagonal matrix. Therefore, each

complex gain element is estimated as
D −1 D
    
∗
ĝp k = [Wi,p ]k ,k X∗i,p k ,k . (34)
i=1 i=1
Finally, the proposed algorithm for the structured Jones
matrices case regarding 3DC calibration regime is given in
2) Estimation of αi : We first need to estimate ϕi,p (we recall Algorithm 2.
that ϕi,p = ηi up + ζi vp ). This is done as follows
ϕ̂i,p = argmin κ̃(ϕi,p ) (35) VI. NUMERICAL SIMULATIONS
ϕi , p
In this part, we first aim to assess the statistical performance
where κ̃(ϕi,p ) = ||Ĵi,p − Gp Hi,p Zi,p (ϕi,p )Fi ||2F . Taking the of Algorithm 1, when the noise model matches our noise as-
derivative of κ̃(ϕi,p ) w.r.t. ϕi,p and setting the result to zero, we sumption, i.e., a SIRP noise modeling. Afterwards, we intend to
obtain study our proposed algorithm, in a more realistic scenario where
 
Tr j exp−j ϕ̂ i , p Ĵi,p FH H H j ϕ̂ i , p
Gp Hi,p Fi ĴH outliers are present. More specifically, the sky is composed of D
i Hi,p Gp − j exp i,p
known bright calibrator sources but also D
weak non-calibrator
=0 (36) sources which are absorbed in the noise component and act as
outliers. Under such assumption, we compare our scheme with
which leads to the recently introduced robust calibration approach based on
 
  Tr Mi,p Student’s t [35] and with the traditional Gaussian cases [32]. Fi-
exp 2j ϕ̂i,p =   (37) nally, we apply Algorithm 2 to the introduced 3DC calibration
Tr MH
i,p
regime where Jones matrices are structured.

where Mi,p = Ĵi,p FH H H

i Hi,p Gp . A. Numerical Results Under SIRP Noise
In the case of a compact array, we can write for the i-th source
The unknown parameters to estimate, θ given in (19), corre-
ϕ̂Ti = α̂Ti Λ (38) spond to the real and imaginary parts of the entries of all Jones
  matrices. The additive noise in (4) is assumed to follow a SIRP
T u1 , . . . , uM
where ϕ̂i = [ϕ̂i,1 , . . . , ϕ̂i,M ] and Λ = . as in (5) and the number of Monte Carlo runs is set to 100.
v1 , . . . , vM
In order to evaluate the estimation performance of the re-
Therefore, estimation of the directional shifts due to propagation
laxed concentrated ML based calibration algorithm, we fix the
in the ionosphere is given by
 M  noise distribution, e.g., as a Student’s t and make use of the

T H p=1 vp
2
− Mp=1 up vp
Cramér-Rao bound (CRB) [65]. To do so, each random texture
ϕ̂i Λ  M component is supposed to follow an inverse gamma distribution
T
− M p=1 vp up p=1 up
2
α̂i = M  M  M
. (39) [66]. As an example,
p=1 up p=1 vp − ( p=1 up vp )
2 2 2

τpq ∼ IG(ν/2, ν/2), (41)

3) Estimation of ϑi : We consider the following minimization
problem with ν degrees of freedom [67] and we choose for example
π
M
 [Ω]k ,l = σ 2 0.9|k −l| expj 2 (k −l) .
ϑ̂i = argmin ||Ĵi,p − Gp Hi,p Zi,p Fi (ϑi )||2F . (40) The covariance inequality principle states that, under quite
ϑi p=1 general/weak conditions, the variance satisfies
 2 
We assume a large number of antennas M while the number
of calibrator sources D is relatively reduced, such that observa- MSE([θ̂]k ) = E [θ̂]k − [θ]k  [CRB(θ)]k ,k (42)
tions outnumber unknown parameters. For each source, the 1D
optimization in (40) can be computed in a reasonable compu- where the CRB is given as the inverse of the Fisher information
tational time through a classical data grid search or a Newton matrix (FIM) F. A Slepian-Bangs type formula of the FIM for
type algorithm. SIRP observations is given in [68] which can be adapted to our
5656 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 21, NOVEMBER 1, 2017

Fig. 2. (a) MSE of the real part of the first 32 unknown parameters for a
given SNR, (b) MSE vs. SNR for the real part of a given unknown parame- Fig. 3. h{θ} as function of the h-th iteration, for loop in line 2 (a), and in
ter and the corresponding CRB, for D = 2 bright signal sources and M = 8 line 1 (b) from Algorithm 1.
antennas, leading to 128 real unknown parameters of interest to estimate and
224 measurements.
Second, we now aim to investigate numerically the conver-
case and reads gence properties of our algorithm, which is composed of 3 loops.
  In each of these loops, θ is updated at each iteration. We there-
ν+4 H
∂ ṽpq (θ) −1 ∂ ṽpq (θ) fore consider the following quantity
[F]k ,l =2  Ω . (43)
ν + 5 pq ∂[θ]k ∂[θ]l

h{θ} = || θ h − θ h−1 ||22 (44)

Note that the noise parameters are decoupled from the param-
where h refers to the h-th iteration. In Fig. 3 we present the
eters of interest. Consequently, only the part corresponding to
convergence rate of loops described in Algorithm 1 at line 1
the latter is kept in the FIM expression.
and 2 (the analysis of convergence of the third loop, given in
First, in Fig. 2(a), we plot the mean square error (MSE)
line 5, has the same behavior as the loop in line 2 and thus, is not
of the real part of each unknown parameter, obtained with
reported here). We note that around 5 iterations are required in
Algorithm 1, for a signal-to-noise ratio SNR = 15 dB. We only
loop 2 to attain convergence while approximately 20 iterations
plot the parameters relative to one given source, the behavior
are needed for the algorithm to be stable, in loop 1. Nevertheless,
being the same for any source. We also compare the MSE of one
in simulations, we notice that only 3 to 4 iterations are sufficient
given parameter, as a function of the SNR, to its corresponding
to get close to the CRB.
CRB in Fig. 2(b). This enables to assess the statistical perfor-
mance of Algorithm 1, and we notice that the MSE approaches
B. Numerical Results Using a Realistic Model
the CRB. The small gap between the bound and the algorithm
is explained by the relaxed nature hypothesis used in the de- Let us investigate the robustness of our proposed calibra-
sign of Algorithm 1. Indeed, we have assumed unknown and tion procedure in a realistic situation, and compare it with the
deterministic texture parameters when we derived the estimates state-of-the-art. To do so, we consider D calibrator sources,
using Algorithm 1, but in the data model, these parameters are D
weak outlier sources and Gaussian background noise in our
generated as random variables following inverse gamma distri- data model. The real parameters of interest to estimate still cor-
bution and the CRB was derived using the prior of the pdf of respond to the real and imaginary entries of the Jones matrices
the texture. associated to the calibrator sources paths. In the following, we
OLLIER et al.: ROBUST CALIBRATION OF RADIO INTERFEROMETERS IN NON-GAUSSIAN ENVIRONMENT 5657

Fig. 5. (a) MSE of the real part of the 16 complex gains for a given SNR,
Fig. 4. (a) MSE of the real part of the 64 unknown parameters for a given SNR, (b) MSE of ζ1 vs. SNR, for D = 2, M = 8, and D
= 4, leading to 38 real
(b) MSE vs. SNR for the real part of a given unknown parameter, for D = 2, parameters of interest to estimate and 224 measurements.
M = 8 and D
= 8, leading to 128 real parameters of interest to estimate and
224 measurements.
source offset ζ1 , cf. Fig 5(b), due to lack of space, the behavior
compare Algorithm 1 with i) the calibration approach exposed in being the same for the other parameters, i.e., ϑ1 , ϑ2 , η2 , η1
[35] which assumes a Student’s t noise modeling using the so- and ζ2 . In the case of structured Jones matrices, adapted to the
called expectation-conditional maximization either algorithm 3DC calibration regime, and in the presence of outliers, we
[69], [70] and ii) the traditional calibration scheme based on still notice the better performances of Algorithm 1, compared
zero-mean white Gaussian noise modeling using a least squares to the state-of-the-art. This is expected since better estimation
approach [32]. The comparative results are plotted in Fig. 4, of the Jones entries leads to better estimation of the physical
for the same computation times. We notice better estimation parameters describing the structured Jones matrices.
performance with Algorithm 1, since we did not specify any
particular noise distribution in our procedure and a SIRP in- VII. CONCLUSION
cludes many different types of distributions. Furthermore, no
assumption was made about independent entries in the noise In this paper, we proposed a robust calibration technique
vector, thus ensuring more flexibility and robustness. where perturbation effects are modeled thanks to Jones ma-
trices. To deal with the presence of outliers in our data, the
C. Structured Case introduced ML estimation method is based on SIRP noise mod-
eling, leading to a relaxed concentrated ML based calibration
In order to compare our proposed global algorithm algorithm. Numerical simulations show that the proposed algo-
(Algorithm 1 followed by Algorithm 2) with the approach based rithm is more robust to the presence of outliers in comparison
on Student’s t [35] and the Gaussian case [32] which were both with the state-of-the-art, for both non-structured and structured
introduced in the non-structured case, we apply Algorithm 2 on Jones matrices, with a reasonable computational complexity.
the output of these two latter algorithms.
Each Jones matrix is generated according to (7) with g =
APPENDIX A
[g1T , . . . , gM
T T
] . In all exposed simulations, we consider similar
computation times for the three presented patterns. We show We describe here the corresponding expressions of (12)
the results only for the complex gains, cf. Fig. 5(a), and the and (16) when we assume a different Ωpq for p < q, p, q ∈
5658 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 21, NOVEMBER 1, 2017

{1, . . . , M }. In this case, the log-likelihood function is written where wi,p = [wi,p(p+1) T T
, . . . , wi,pM ]T , ui,p (θ i,p ) = [uTi,p(p+1)
as (θ i,p ), . . . , uTi,pM (θ i,p )]T and Ai,p = bdiag{βi τp(p+1) Ω, . . . ,
log f (x|θ, τ , Ω12 , Ω13 , . . . , Ω(M −1)M ) = −4B log π βi τpM Ω}−1 .
∗T ∗T
   1 Furthermore, we have w̃i,p = [wi,1p , . . . , wi,(p−1)p ]T ,
−1
−4 log τpq − log |Ωpq | − aH
pq (θ)Ωpq apq (θ).
T T
ũi,p (θ i,p ) = [u∗i,1p (θ i,p ), . . . , u∗i,(p−1)p (θ i,p )]T and Ãi,p =
pq pq pq
τ pq
bdiag{βi τ1p Ω∗ , . . . , βi τ(p−1)p Ω∗ }−1 .
(45)
We make use of (27) in what follows
For each antenna pair, the texture estimate is given by ⎡ ⎤ ⎡ ⎤
ui,p(p+1) (θ i,p ) Σi,p+1 θ i,p
1 H ⎢ ⎥ ⎢ ⎥
a (θ)Ω−1
τ̂pq = pq apq (θ) (46) ui,p (θ i,p ) = ⎣ ..
⎦=⎣
..
⎦ = Σi θ i,p
4 pq . .
while the speckle covariance estimate reads ui,pM (θ i,p ) Σi,M θ i,p
(51)
h+1 apq (θ)aHpq (θ) where Σi = [ΣTi,p+1 , · · · , ΣTi,M ]T . Likewise, we use (28) in
Ω̂pq =4  h −1 . (47)
aH (θ) Ω̂ a (θ) ⎡ ⎤ ⎡ ⎤
pq pq pq u∗i,1p (θ i,p ) Υ∗i,1 θ i,p
⎢ .. ⎥ ⎢ .. ⎥
The remainder of the algorithm is straightforwardly obtained ũi,p (θ i,p ) = ⎣ . ⎦=⎣ . ⎦ = Υi θ i,p
using (47). u∗i,(p−1)p (θ i,p ) Υ∗i,p−1 θ i,p
(52)
T T
APPENDIX B in which Υi = [Υ∗i,1 , · · · , Υ∗i,p−1 ]T . Inserting (51) and (52)
We present here the steps to obtain (29). Firstly, for sake into (50) and taking the derivative w.r.t. θ i,p leads to the expres-
of clarity, let us denote ci = [ci 1 , ci 2 , ci 3 , ci 4 ]T to refer to the sions in (29), using the fact that Ai,p and Ãi,p are Hermitian.
four entries of the vectorization of source coherency matrix Ci .
Likewise, for the i-th source, we write REFERENCES
 
p pi 2
Ji,p (θ i,p ) = i 1 [1] A. R. Thompson, J. M. Moran, and G. W. Swenson Jr., Interferometry and
pi 3 pi 4 Synthesis in Radio Astronomy, 2nd ed. Hoboken, NJ, USA: Wiley, 2001.
[2] M. P. Van Haarlem et al., “LOFAR: The LOw-Frequency ARray,” Astron.
for the p-th antenna and Astrophys., vol. 556, no. A2, 2013, 53 pages.
  [3] P. E. Dewdney, P. J. Hall, R. T. Schilizzi, and T. J. L. Lazio, “The Square
q qi 2 Kilometre Array,” Proc. IEEE, vol. 97, no. 8, pp. 1482–1496, Aug. 2009.
Ji,q (θ i,q ) = i 1
qi 3 qi 4 [4] S. J. Wijnholds and A.-J. van der Veen, “Fundamental imaging limits of
radio telescope arrays,” IEEE J. Sel. Topics Signal Process., vol. 2, no. 5,
for the q-th antenna, i.e., θ i,p = [pi 1 , pi 2 , pi 3 , pi 4 ]T and θ i,q = pp. 613–623, Oct. 2008.
[5] U. Rau, S. Bhatnagar, M. A. Voronkov, and T. J. Cornwell, “Advances in
[qi 1 , qi 2 , qi 3 , qi 4 ]T . Using these latter notation, we obtain (27) calibration and imaging techniques in radio interferometry,” Proc. IEEE,
where vol. 97, no. 8, pp. 1472–1481, Aug. 2009.
⎡ ⎤ [6] M. de Vos, A. W. Gunst, and R. Nijboer, “The LOFAR telescope: System
αi,q βi,q 0 0 architecture and signal processing,” Proc. IEEE, vol. 97, no. 8, pp. 1431–
⎢ 0 αi,q βi,q ⎥ 1437, Aug. 2009.
⎢ 0 ⎥
Σi,q = ⎢ ⎥ (48) [7] A.-J. van der Veen and S. J. Wijnholds, “Signal processing tools for radio
⎣ γi,q ρi,q 0 0 ⎦ astronomy,” in Handbook of Signal Processing Systems. New York, NY,
USA: Springer-Verlag, 2013, pp. 421–463.
0 0 γi,q ρi,q [8] D. A. Mitchell et al., “Real-time calibration of the Murchison Widefield
in which αi,q = qi∗1 ci 1 + qi∗2 ci 3 , βi,q = qi∗1 ci 2 + qi∗2 ci 4 , γi,q = Array,” IEEE J. Sel. Topics Signal Process., vol. 2, no. 5, pp. 707–717,
Oct. 2008.
qi∗3 ci 1 + qi∗4 ci 3 and ρi,q = qi∗3 ci 2 + qi∗4 ci 4 . [9] S. J. Wijnholds, A.-J. van der Veen, F. De Stefani, E. La Rosa, and A.
We also obtain (28) where Farina, “Signal processing challenges for radio astronomical arrays,” in
⎡ ⎤ Proc. Int. Conf. Acoust., Speech, Signal Process., Florence, Italy, 2014,
λi,q μi,q 0 0 pp. 5382–5386.
⎢ν 0 ⎥
[10] S. J. Wijnholds, S. van der Tol, R. Nijboer, and A.-J. van der Veen,
⎢ i,q ξi,q 0 ⎥
Υi,q = ⎢ ⎥ (49) “Calibration challenges for future radio telescopes,” IEEE Signal Process.
⎣ 0 0 λi,q μi,q ⎦ Mag., vol. 27, no. 1, pp. 30–42, Jan. 2010.
[11] D. R. Fuhrmann, “Estimation of sensor gain and phase,” IEEE Trans.
0 0 νi,q ξi,q Signal Process., vol. 42, no. 1, pp. 77–87, Jan. 1994.
[12] B. C. Ng and C. M. S. See, “Sensor-array calibration using a maximum-
in which λi,q = qi 1 ci 1 + qi 2 ci 2 , μi,q = qi 1 ci 3 + qi 2 ci 4 , νi,q = likelihood approach,” IEEE Trans. Antennas Propag., vol. 44, no. 6,
qi 3 ci 1 + qi 4 ci 2 and ξi,q = qi 3 ci 3 + qi 4 ci 4 . pp. 827–835, Jun. 1996.
Finally, the cost function in (26) can be written as [13] B. P. Ng, J. P. Lie, M. H. Er, and A. Feng, “A practical simple geometry
and gain/phase calibration technique for antenna array processing,” IEEE
 H   Trans. Antennas Propag., vol. 57, no. 7, pp. 1963–1972, Jul. 2009.
φi (θ i,p ) = wi,p − ui,p (θ i,p ) Ai,p wi,p − ui,p (θ i,p ) [14] J. T. H. Lo and S. L. Marple, “Eigenstructure methods for array sensor
localization,” in Proc. Int. Conf. Acoust., Speech, Signal Process., Dallas,
 H   TX, USA, 1987, pp. 2260–2263.
+ w̃i,p − ũi,p(θ i,p ) Ãi,p w̃i,p − ũi,p (θ i,p ) + Constant [15] B. C. Ng and W. Ser, “Array shape calibration using sources in known loca-
tions,” in Proc. ICCS/ISITA ‘Commun. Move,’ Singapore, 1992, pp. 836–
(50) 840.
OLLIER et al.: ROBUST CALIBRATION OF RADIO INTERFEROMETERS IN NON-GAUSSIAN ENVIRONMENT 5659

[16] B. C. Ng and A. Nehorai, “Active array sensor localization,” Signal Pro- [42] J. P. Hamaker, J. D. Bregman, and R. J. Sault, “Understanding radio
cess., vol. 44, no. 3, pp. 309–327, 1995. polarimetry. I. Mathematical foundations,” Astron. Astrophys. Suppl. Ser.,
[17] Y. Rockah and P. Schultheiss, “Array shape calibration using sources vol. 117, no. 1, pp. 137–147, 1996.
in unknown locations—Part I: Far-field sources,” IEEE Trans. Acoust., [43] O. M. Smirnov, “Revisiting the radio interferometer measurement equa-
Speech, Signal Process., vol. 35, no. 3, pp. 286–299, Mar. 1987. tion. I. A full-sky Jones formalism,” Astron. Astrophys., vol. 527, no. A106,
[18] A. J. Weiss and B. Friedlander, “Array shape calibration using sources 11 pages, 2011.
in unknown locations-a maximum likelihood approach,” IEEE Trans. [44] J. E. Noordam, “The measurement equation of a generic radio telescope,
Acoust., Speech, Signal Process., vol. 37, no. 12, pp. 1958–1966, AIPS++ implementation note nr 185,” ASTRON, Dwingeloo, The Nether-
Dec. 1989. lands, Tech. Rep., 1996.
[19] B. Friedlander and A. J. Weiss, “Direction finding in the presence of mutual [45] J. E. Noordam and O. M. Smirnov, “The MeqTrees software system and
coupling,” IEEE Trans. Antennas Propag., vol. 39, no. 3, pp. 273–284, its use for third-generation calibration of radio interferometers,” Astron.
Mar. 1991. Astrophys., vol. 524, no. A61, 16 pages, 2010.
[20] M. P. Wylie, S. Roy, and H. Messer, “Joint DOA estimation and phase cal- [46] C. Lonsdale, “Calibration approaches,” LFD memo 015, Tech. Rep.,
ibration of linear equispaced (LES) arrays,” IEEE Trans. Signal Process., 2004.
vol. 42, no. 12, pp. 3449–3459, Dec. 1994. [47] L. M. Ker, “Radio AGN evolution with low frequency radio surveys,”
[21] L. Qiong, G. Long, and Y. Zhongfu, “An overview of self-calibration in Ph.D. dissertation, Univ. Edinburgh, Edinburgh, U.K., 2012.
sensor array processing,” in Proc. 6th Int. Symp. Antennas, Propag. EM [48] C. D. Nunhokee, “Link between ghost artefacts, source suppression and
Theory, Beijing, China, 2003, pp. 279–282. incomplete calibration sky models,” Ph.D. dissertation, Dept. Phys. Elec-
[22] B. P. Flanagan and K. L. Bell, “Array self-calibration with large sensor tron., Rhodes Univ., Rhodes, South Africa, 2015.
position errors,” Signal Process., vol. 81, no. 10, pp. 2201–2214, 2001. [49] X. Zhang, M. N. El Korso, and M. Pesavento, “MIMO radar target local-
[23] J. Baars, R. Genzel, I. Pauliny-Toth, and A. Witzel, “The absolute spectrum ization and performance evaluation under SIRP clutter,” Signal Process.
of CAS A-an accurate flux density scale and a set of secondary calibrators,” J., vol. 130, pp. 217–232, 2017.
Astron. Astrophys., vol. 61, no. 1, pp. 99–106, 1977. [50] S. Yatawatta, “Reduced ambiguity calibration for LOFAR,” Exp. Astron-
[24] A.-J. Boonstra and A.-J. van der Veen, “Gain calibration methods for radio omy, vol. 34, no. 1, pp. 89–103, 2012.
telescope arrays,” IEEE Trans. Signal Process., vol. 51, no. 1, pp. 25–38, [51] S. Kazemi, P. Hurley, O. Öçal, and G. Cherubini, “Blind calibration for
Jan. 2003. radio interferometry using convex optimization,” in Proc. 3rd Int. Work-
[25] S. J. Wijnholds and A.-J. van der Veen, “Multisource self-calibration for shop Compressed Sens. Theory Appl. Radar, Sonar Remote Sens., Pisa,
sensor arrays,” IEEE Trans. Signal Process., vol. 57, no. 9, pp. 3512–3522, Italy, 2015, pp. 164–168.
Sep. 2009. [52] S. van der Tol, B. D. Jeffs, and A.-J. van der Veen, “Self-calibration for the
[26] S. J. Wijnholds, “Fish-eye observing with phased array radio telescopes,” LOFAR radio astronomical array,” IEEE Trans. Signal Process., vol. 55,
Ph.D. dissertation, Dept. Microelectron. Comput. Eng., Technische Uni- no. 9, pp. 4497–4510, Sep. 2007.
versiteit Delft, Delft, The Netherlands, 2010. [53] E. Conte, A. De Maio, and G. Ricci, “Recursive estimation of the co-
[27] S. Salvini and S. J. Wijnholds, “Fast gain calibration in radio astronomy variance matrix of a compound-Gaussian process and its application to
using alternating direction implicit methods: Analysis and applications,” adaptive CFAR detection,” IEEE Trans. Signal Process., vol. 50, no. 8,
Astron. Astrophys., vol. 571, no. A97, 14 pages, 2014. pp. 1908–1915, Aug. 2002.
[28] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from [54] F. Pascal, Y. Chitour, J.-P. Ovarlez, P. Forster, and P. Larzabal, “Covariance
incomplete data via the EM algorithm,” J. Roy. Statist. Soc. B (Methodol.), structure maximum-likelihood estimates in compound Gaussian noise:
vol. 39, no. 1, pp. 1–38, 1977. Existence and algorithm analysis,” IEEE Trans. Signal Process., vol. 56,
[29] T. K. Moon, “The expectation-maximization algorithm,” IEEE Signal no. 1, pp. 34–48, Jan. 2008.
Process. Mag., vol. 13, no. 6, pp. 47–60, Nov. 1996. [55] X. Zhang, M. N. El Korso, and M. Pesavento, “Maximum likelihood and
[30] G. J. McLachlan and T. Krishnan, The EM Algorithm and Extensions, 2nd maximum a posteriori direction-of-arrival estimation in the presence of
ed. Hoboken, NJ, USA: Wiley, 2008. SIRP noise,” in Proc. Int. Conf. Acoust., Speech, Signal Process., Shang-
[31] J. A. Fessler and A. O. Hero, “Space-alternating generalized expectation- hai, China, 2016, pp. 3081–3085.
maximization algorithm,” IEEE Trans. Signal Process., vol. 42, no. 10, [56] V. Ollier, M. N. El Korso, R. Boyer, P. Larzabal, and M. Pesavento,
pp. 2664–2677, Oct. 1994. “Relaxed concentrated MLE for robust calibration of radio interferome-
[32] S. Yatawatta, S. Zaroubi, G. de Bruyn, L. Koopmans, and J. Noordam, ters,” in Proc. 24th Eur. Signal Process. Conf., Budapest, Hungary, 2016,
“Radio interferometric calibration using the SAGE algorithm,” in Proc. pp. 280–284.
13th Digit. Signal Process. Workshop/5th Signal Process. Edu. Workshop, [57] A. Hjorungnes and D. Gesbert, “Complex-valued matrix differentiation:
Marco Island, FL, USA, 2009, pp. 150–155. Techniques and key results,” IEEE Trans. Signal Process., vol. 55, no. 6,
[33] J. Raza, A.-J. Boonstra, and A.-J. Van der Veen, “Spatial filtering of RF pp. 2740–2746, Jun. 2007.
interference in radio astronomy,” IEEE Signal Process. Lett., vol. 9, no. 2, [58] M. Feder and E. Weinstein, “Parameter estimation of superimposed signals
pp. 64–67, Feb. 2002. using the EM algorithm,” IEEE Trans. Acoust., Speech, Signal Process.,
[34] S. Van Der Tol and A.-J. Van Der Veen, “Performance analysis of spa- vol. 36, no. 4, pp. 477–489, Apr. 1988.
tial filtering of RF interference in radio astronomy,” IEEE Trans. Signal [59] K. Madsen, H. B. Nielsen, and O. Tingleff, “Methods for non-linear least
Process., vol. 53, no. 3, pp. 896–910, Mar. 2005. squares problems,” Dept. Informat. Math. Model., Tech. Univ. Denmark,
[35] S. Yatawatta and S. Kazemi, “Robust radio interferometric calibration,” Kgs. Lyngby, Denmark, Tech. Rep., 2004.
in Proc. Int. Conf. Acoust., Speech, Signal Process., Florence, Italy, 2014, [60] J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. New York,
pp. 5429–5433. NY, USA: Springer-Verlag, 2006.
[36] A.-J. Boonstra, “Radio frequency interference mitigation in radio astron- [61] H. Gavin, “The Levenberg-Marquardt method for nonlinear least squares
omy,” Ph.D. dissertation, Elect. Eng., Math. and Comput. Sci. Faculty, curve-fitting problems,” Duke Univ., Durham, NC, USA, 2011.
Technische Univ. Delft, Delft, The Netherlands, 2005. [62] T. W. Anderson, An Introduction to Multivariate Statistical Analysis,
[37] E. Jay, “Détection en environnement non gaussien,” Ph.D. dissertation, vol. 2. New York, NY, USA: Wiley, 1958.
ETIS—Equipe traitement des images et du signal, Univ. de Cergy Pon- [63] D. P. Bertsekas, Nonlinear Programming. Belmont, MA, USA: Athena
toise, Cergy, France, 2002. Scientific, 1999.
[38] K. Yao, “Spherically invariant random processes: Theory and ap- [64] S. Vorobyov, A. B. Gershman, and K. M. Wong, “Maximum likelihood
plications,” in Communications, Information and Network Security, direction-of-arrival estimation in unknown noise fields using sparse sen-
V. Bhargava et al., Eds. Dordrecht, The Netherlands: Kluwer, 2003, sor arrays,” IEEE Trans. Signal Process., vol. 53, no. 1, pp. 34–43,
pp. 315–331. Jan. 2005.
[39] J. Wang, A. Dogandzic, and A. Nehorai, “Maximum likelihood estimation [65] P. Stoica and R. L. Moses, Spectral Analysis of Signals. Upper Saddle
of compound-Gaussian clutter and target parameters,” IEEE Trans. Signal River, NJ, USA: Prentice-Hall, 2005.
Process., vol. 54, no. 10, pp. 3884–3898, Oct. 2006. [66] A. Balleri, A. Nehorai, and J. Wang, “Maximum likelihood estima-
[40] J. Friedman et al., “Pathwise coordinate optimization,” Ann. Appl. Statist., tion for compound-Gaussian clutter with inverse gamma texture,” IEEE
vol. 1, no. 2, pp. 302–332, 2007. Trans. Aerosp. Electron. Syst., vol. 43, no. 2, pp. 775–780, Apr.
[41] M. Hong, M. Razaviyayn, Z.-Q. Luo, and J.-S. Pang, “A unified algo- 2007.
rithmic framework for block-structured optimization involving big data: [67] E. Ollila, D. E. Tyler, V. Koivunen, and H. Poor, “Complex elliptically
With applications in machine learning and signal processing,” IEEE Signal symmetric distributions survey, new results and applications,” IEEE Trans.
Process. Mag., vol. 33, no. 1, pp. 57–77, Jan. 2016. Signal Process., vol. 60, no. 11, pp. 5597–5625, Nov. 2012.
5660 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 65, NO. 21, NOVEMBER 1, 2017

[68] O. Besson and Y. I. Abramovich, “On the Fisher information matrix for Pascal Larzabal (M’93) was born in Basque coun-
multivariate elliptically contoured distributions,” IEEE Signal Process. try in the south of France in 1962. In 1985, he joined
Lett., vol. 20, no. 11, pp. 1130–1133, Nov. 2013. the Ecole Normale Supérieure of Cachan, Cachan,
[69] C. Liu and D. B. Rubin, “ML estimation of the t distribution using EM and France, from where he received the Agrégation in
its extensions, ECM and ECME,” Statistica Sinica, vol. 5, no. 1, pp. 19–39, electrical engineering in 1988. He received the Ph.D.
1995. degree in 1992 and the Habilitation à Diriger des
[70] S. Li, H. Wang, and T. Chai, “A t-distribution based particle filter for target Recherches in 1998. From 1998 to 2003, he was
tracking,” in Proc. Amer. Control Conf., Minneapolis, MN, USA, 2006, the Director of Electrical Engineering IUP, Univer-
pp. 2191–2196. sity Paris-Sud. From March 2003 to March 2007, he
was at the head of the Electrical Engineering De-
partment, Institut Universitaire de Technologie of
Virginie Ollier was born in Beaumont, France, in Cachan. He is currently a Professor of electrical engineering at the Univer-
1991. She received the Master Research degree in sig- sity of Paris-Sud XI. He teaches electronic, signal processing, control, and
nal and image processing from Ecole Centrale Mar- mathematics. Since January 2007, he has been the Director of the Laboratory
seille, Marseille, France, in 2015. She is currently Systèmes et Applications des Technologies de l’Information et de l’Energie,
working toward the Ph.D. degree in signal process- Paris. His research concerns estimation in array processing and spectral analy-
ing at the Ecole Normale Supérieure Paris-Saclay, sis for wavefront identification, radars, communications, tomography, and med-
Cachan, France, and is a member of Systèmes et Ap- ical imaging. His recent works concern geographical positioning and radio
plications des Technologies de l’Information et de astronomy.
l’Energie Laboratory and the Laboratoire des Signaux
et Systèmes. Her research interests include estimation
theory and statistical signal processing, especially for
applications in robust array signal processing.

Mohammed Nabil El Korso was born in Oran,

Algeria. He received the M.Sc. degree in electrical
engineering from the National Polytechnic School,
El-Harrach, Algeria, in 2007, and the Master Re-
search degree in signal and image processing and
the Ph.D. degree, both from Paris-Sud XI University,
France, in 2008 and 2011, respectively. From 2011
to 2012, he was a Research Scientist in the Com-
munication Systems Group at Technische Universität
Darmstadt, Germany. He was an Assistant Professor
at Ecole Normale Supérieure de Cachan from 2012
to 2013. He is currently an Assistant Professor at the University of Paris Ouest
Nanterre la Defense and a member of LEME (EA4416) Laboratory. His re-
search interests include statistical signal processing and estimation/detection
theory with applications to array signal processing.

Rémy Boyer (SM’16) received the M.Sc. and Ph.D.

degrees from the Ecole Nationale Supérieure des
Telecommunications (ENST-Paris or Télécom Paris-
Tech), Paris, France, in 1999 and 2002, respectively,
in statistical signal processing. From 2002 to 2003,
he was a Postdoctoral Fellow for six months at Sher-
brooke University (Canada). Since 2003, he has been Marius Pesavento (M’00) received the Dipl.Ing. de-
an Associate Professor in the Laboratory of Signals gree in 1999 and the Dr. Ing. degree in electrical engi-
and Systems (L2S), University of Paris-Sud. From neering in 2005, both from Ruhr-University Bochum,
2011 to 2012, he was a visiting researcher at Systèmes Bochum, Germany, and the M.Eng. degree from Mc-
et Applications des Technologies de l’Information et Master University, Hamilton, ON, Canada, in 2000.
de l’Energie Laboratory (Ecole Normale Supérieure de Cachan) and at the Between 2005 and 2009, he held research positions in
University of Aalborg (Denmark). His research interests include compressive two start-up companies. In 2010, he became an As-
sampling, array signal processing, Bayesian performance bounds for parameter sistant Professor for robust signal processing and a
estimation and detection, security in mobile networks, as well as numerical linear Full Professor for communication systems in 2013 in
and multilinear algebra. He received the “Habilitation à Diriger des Recherches” the Department of Electrical Engineering and Infor-
from the University of Paris-Sud in December 2012. mation Technology, Technical University Darmstadt,
Rémy Boyer is the author or co-author of over 120 publications. He was Darmstadt, Germany. His research interests include robust signal processing and
awarded the prize for excellence in scientific research (“prime d’excellence adaptive beamforming, high-resolution sensor array processing, multiantenna
scientifique” in french) for 12 years. He is an elected member of the CCSU and multiuser communication systems, distributed, sparse, and mixed-integer
(“Commission Consultative de Spécialistes de l’Université”), an elected mem- optimization techniques for signal processing and communications, statistical
ber of the L2S laboratory board, a nominated member of the CNU (“Conseil signal processing, spectral analysis, and parameter estimation. He received the
National des Universités”) in section 61 and an expert for the ANR (“Agence 2003 ITG/VDE Best Paper Award, the 2005 Young Author Best Paper Award
Nationale de la Recherche”). Rémy Boyer is a member of the EURASIP Spe- of the IEEE TRANSACTIONS ON SIGNAL PROCESSING, and the 2010 Best Paper
cial Area Team SPMuS (Signal Processing for Multisensor Systems), an affiliate Award of the CrownCOM Conference. He is a Member of the Editorial Board
member of the IEEE Sensor Array and Multichannel (SAM) Technical Com- of the EURASIP Signal Processing Journal, and served as an Associate Editor
mittee. Rémy Boyer is also the elected vice-president/secretary of the French of the IEEE TRANSACTIONS ON SIGNAL PROCESSING during 2012–2016. He is
chapter of the IEEE Signal Processing Society, the scientific director of the 12-th a Member of the Sensor Array and Multichannel Technical Committee of the
French summer school in Signal and Image Processing and an elected member IEEE Signal Processing Society, and the Special Area Teams “Signal Processing
of the board of the national research group (“GDR”) ISIS (Information, Signal, for Communications and Networking” and “Signal Processing for Multisensor
Image and viSion). Systems” of the EURASIP.