Je Pense

Gwenaël, Hélène, Sophie, Béatrice, Aurore, Dorothée, Fabien Jean-Pierre Dominique

Muriel, Yohan, Mehdi, Alexiane, Sakina mais aussi David

Alex

Lucien

Yassine Gordana

Arnaud

Jean-Phi Rayen

Chengfang

Ammar, Eugénie Eric C Bruno

Abigael Martin

Alice Cyrile

Abdellatif

Maxime, Julien Guillaume

autres... J'aimerais Tant D'

Nhung Philippe Yannick

Marie Et

Virginie Ollier

Keywords: Parametric estimation, calibration, robustness, distributed optimization, radio astronomy, array processing. x Estimation paramétrique, calibration, robustesse, optimisation distribuée, radioastronomie, traitement d'antennes. xii

Une thèse, c'est la chance de pouvoir travailler pendant trois ans sur un sujet qui nous tient à coeur et acquérir ainsi une expertise technique qui constitue un atout essentiel d'un point de vue professionnel. Mais c'est aussi des échanges et des rencontres avec une multitude de personnes qui font de la thèse une véritable aventure humaine. Je tiens tout d'abord à remercier l'ensemble des membres du jury pour l'intérêt qu'ils ont porté à mon travail et leurs remarques qui ont permis de l'enrichir. Je remercie David Brie et Jean-François Giovannelli d'avoir accepté d'être rapporteurs sur ma thèse et qui m'ont apporté des commentaires très constructifs et utiles. Merci à Sylvie Marcos et Jean-Luc Starck d'avoir accepté d'évaluer mes travaux lors du suivi à mi-parcours, avec Jean-François Giovannelli également, et qui étaient aussi présents lors de ma soutenance en tant qu'examinateurs, avec des remarques toujours pertinentes. Je remercie Christian Jutten de m'avoir fait l'honneur de présider ce comité et d'apporter son avis sur ce travail. Selon moi, la qualité d'une thèse se mesure essentiellement par l'encadrement. Et de ce point de vue là, j'ai été très chanceuse. Merci Pascal pour tous les bons conseils que tu as pu me donner, pour ton enthousiasme sans faille, ton expérience sur un très grand nombre de sujets et ton soutien qui m'a permis de garder conance tout le long de la thèse. Merci aussi de m'avoir fait conance et de m'avoir donné l'opportunité de faire cette thèse après le stage, en me laissant aussi la chance de voyager vers des destinations incroyables lors de conférences ou écoles thématiques. Merci pour tout.

Merci Rémy pour ton aide très précieuse, ta très forte capacité de travail et ta très grande disponibilité. J'ai aussi eu la chance de pouvoir eectuer mes tout premiers enseignements à tes côtés et ce fut une expérience très agréable. Tu es toujours de très bon conseil avec une très grande rigueur scientique et j'espère vraiment que tu te plairas à Lille ! Merci Nabil. Tu as toujours été très présent, que ce soit pendant le stage ou la thèse, très disponible, investi et à l'écoute. Tu as pu m'aider au jour le jour, répondre à toutes les questions que je pouvais avoir, me donner de nombreux conseils et m'accorder beaucoup de temps avec une patience et un soutien indéfectibles. Merci pour cela, j'ai vraiment apprécié pouvoir travailler à tes côtés.

Grâce à vous, cette thèse fut pour moi une expérience inoubliable et j'espère que nous pourrons toujours garder contact et nous revoir dans les années à venir.

Une thèse, c'est aussi un environnement de recherche privilégié au sein de laboratoires : le Satie et le L2S dans mon cas. Je tiens à remercier l'ensemble des personnes que j'ai pu côtoyer de près ou de loin et qui rendent l'atmosphère de travail toujours très chaleureuse et agréable, de par leur sympathie et leur bonne humeur.

Résumé

Une réelle mutation technologique s'opère dans la conception des instruments d'observation pour la radioastronomie. En eet, les nouveaux réseaux phasés à balayage électronique viennent supplanter les traditionnelles antennes paraboliques à balayage mécanique et s'appuyent sur l'interférométrie pour obtenir des images de meilleure résolution. Les signaux radio mesurés par chaque antenne sont alors corrélés entre eux. Grâce à la synthèse d'ouverture, la résolution angulaire résultante est équivalente à celle obtenue avec un unique large télescope dont le diamètre équivaudrait à la taille du réseau. Ces nouveaux réseaux d'antennes sont particulièrement exibles puisqu'ils sont orientés électroniquement selon diérentes directions possibles, sans aucun mouvement mécanique des antennes. Parmi la nouvelle génération d'instruments très sensibles, nous pouvons mentionner le LOFAR (Low Frequency Array) qui regroupe près de 50 000 antennes à travers l'Europe mais aussi le futur radiotélescope géant SKA (Square Kilometre Array) dont les stations seront réparties à terme en Australie et en Afrique du Sud et fourniront de nouvelles études essentielles sur notre galaxie. Néanmoins, la grande surface collectrice du réseau, le très large champ de vision des antennes et la forte dynamique des images nécessitent de relever un certain nombre de dés scientiques en termes de calibration, réduction des données et reconstruction d'images.

En radioastronomie, les signaux d'intérêt sont perturbés par de nombreux eets environnementaux et instrumentaux, nécessitant la mise en oeuvre de techniques algorithmiques pour les traiter et pouvoir ainsi reconstruire in ne des images parfaitement nettes de l'espace. Cette étape de correction des perturbations se nomme la calibration et repose généralement sur une méthode des moindres carrés, en s'appuyant sur des sources de référence dont la position et l'intensité sont connues. Cette méthode classique se révèle être théoriquement asymptotiquement ecace dans le cas d'un bruit Gaussien et fonctionne généralement de manière itérative pour une seule fréquence considérée. Cependant, en pratique, l'hypothèse du bruit Gaussien n'est pas toujours valide car de multiples sources inconnues de faible intensité sont visibles dans le champ de vision et des interférences radioélectriques perturbent les données pour nalement modier la nature typiquement gaussienne du bruit. En outre, réaliser une calibration fréquence par fréquence, de façon complètement indépendante, n'est pas la manière la plus optimale de procéder. Le but de ce travail est donc de développer des algorithmes de correction dans le traitement des signaux radio qui soient robustes à la présence d'éventuelles valeurs aberrantes ou sources d'interférences, et qui soient adaptés au contexte multi-fréquentiel.

Par conséquent, nous nous appuyons sur une modélisation plus générale que la loi gaussienne, appelée processus Gaussien composé, qui inclut un grand nombre de distributions diérentes telles que la loi de Student ou la loi de Laplace. L'algorithme robuste que nous proposons est itératif et basé sur l'estimation au sens du maximum de vraisemblance. Pour réduire les temps de calcul, nous avons recours aux algorithmes EM (Expectation Maximization) et BCD (Block Coordinate Descent) qui, pour un choix judicieux de paramétrisation, fournissent des expressions analytiques pour les paramètres d'intérêt. En accord avec le scénario multi-fréquentiel sous étude, xi nous exploitons la variation spectrale des perturbations en utilisant des méthodologies telles que l'optimisation distribuée sous contraintes et le traitement parallèle des données, notamment avec l'algorithme ADMM (Alternating Direction Method of Multipliers).

3.2

Mean dierence between true and estimated positions for the D background sources. • ADMM: Alternating Direction Method of Multipliers.

• AIPS: Astronomical Image Processing System.

• ASKAP: Australian Square Kilometre Array Pathnder.

• ASTRON: Netherlands Institute for Radio Astronomy.

• BBS: Black-Board Self-cal.

• BCD: Block Coordinate Descent.

• BFGS: Broyden-Fletcher-Goldfarb-Shanno.

• CASA: Common Astronomy Software Applications.

• CasA: Cassiopeia A.

• CG: Compound-Gaussian.

• COLR: Column Ratio.

• CRB: Cramer-Rao Bound.

• CS: Cotton-Schwab.

• CygA: Cygnus A.

• DEC: Declination.

• DOA: Direction Of Arrival.

• ECME: Expectation-Conditional Maximization Either.

• EM: Expectation Maximization.

• FBC: Field-Based Calibration.

• FIM: Fisher Information Matrix.

• GNLS: Generalized Nonlinear Least Squares.

• GPU: Graphics Processing Unit.

• IG-CG: Inverse-Gaussian Compound-Gaussian.

• i.i.d.: independent and identically distributed.

• IMAPE: Iterative MAP Estimator.

xix

• IRLS: Iterative Reweighted Least Squares.

• LM: Levenberg-Marquardt.

• LOFAR: Low Frequency Array.

• LOGLS: Logarithmic Least Squares.

• LS: Least-Squares.

• LS-MVI: Least Squares Minimum Variance Imaging.

• MAP: Maximum A Posteriori.

• MeerKAT: Karoo Array Telescope.

• MEM: Maximum Entropy Method.

• ML: Maximum Likelihood.

• MMSE: Minimum Mean Squared Error.

• MRCA: Multi-frequency Robust Calibration Algorithm.

• MRCA R3 : Multi-frequency Robust Calibration Algorithm for Regime 3.

• MS: Measurement Set.

• MSE: Mean Square Error.

• MWA: Murchison Wideeld Array.

• OS: Ordered Subsets.

• pdf: probability density function.

• PSF: Point Spread Function.

• RA: Right Ascension.

• RCMLE: Relaxed Concentrated ML Estimator.

• RM: Rotation Measure.

• RMT: Random Matrix Theory.

• RTR: Riemannian Trust-Region.

• SA: Simulated Annealing.

• SAGE: Space Alternating Generalized Expectation Maximization.

• SKA: Square Kilometre Array.

• SNR: Signal-to-Noise Ratio.

• SPAM: Source Peeling and Atmospheric Modeling.

• StEFCal: Statistically EFcient and Fast Calibration.

• SUMSS: Sydney University Molonglo Sky Survey.

• TEC: Total Electron Content.

• VLA: Very Large Array.

• (W)ALS: (Weighted) Alternating Least-Squares.

• WSF: Weighted Subspace Fitting.

• WSRT: Westerbork Synthesis Radio Telescope.

xx Mathematical symbols

• C : complex numbers.

• R : real numbers.

• {•} : real part of a complex number.

• Im {•} : imaginary part of a complex number.

• a : scalar quantity.

• a : column vector.

• A : matrix.

• (•) * : complex conjugate operator.

• (•) T : transpose operator.

• (•) H : Hermitian transpose operator.

• I N : N × N identity matrix.

• | • | : determinant operator.

• || • || 2 : l 2 norm.

• || • || F : Frobenius norm.

• [a] k : k-th entry of vector a.

• [A] i,j : entry in the i-th row and j-th column of matrix A.

• [A] :,k : k-th column vector of matrix A.

• tr {•} : trace operator.

• vec (•) : vectorization operator to stack all columns of a matrix on top of one another.

• N (m, C) : Gaussian probability density function with mean m and covariance matrix C.

• p (x) : probability density function.

• p (x, y) : joint probability density function.

• p ( x| y) : conditional probability density function.

• E{•} : expectation operator.

• ⊗ : Kronecker product between two matrices. For two coordinate vectors, we refer to it as the outer product.

• bdiag{•} : block-diagonal operator.

• diag{•} : conversion of a vector into a diagonal matrix.

• j : complex number whose square equals -1.

xxi xxii

Résumé étendu Introduction

Le travail présenté dans ce document porte sur l'étude des nouveaux radio-interféromètres tels que le LOFAR (Low Frequency Array) [START_REF] Van Haarlem | LOFAR: The LOw-Frequency ARray[END_REF] en Europe ou le futur SKA (Square Kilometre Array) [START_REF] Dewdney | The Square Kilometre Array[END_REF] Dans ce document, nous nous intéressons plus particulièrement à l'étape de calibration qui consiste à corriger les perturbations subies par le signal radio depuis son émission par une source jusqu'à sa réception par le réseau d'antennes.

Parmi les eets perturbateurs, nous pouvons citer

• la présence des gains électroniques complexes des antennes.

• l'anisotropie dans la réponse des antennes (réponse non uniforme, d'où une dépendance selon la direction d'arrivée).

• le retard de propagation dû au passage du signal par la troposphère et l'ionosphère. Ces couches spéciques de l'atmosphère se caractérisent par des phénomènes de diraction et de réfraction, variables dans le temps et l'espace, et aectent notamment la polarisation du signal.

En calibration, la plupart des algorithmes s'appuient sur une méthode de type moindres carrés, en minimisant une fonction coût par rapport à chaque paramètre inconnu, de façon alternée [START_REF]Multisource self-calibration for sensor arrays[END_REF][START_REF] Boonstra | Gain calibration methods for radio telescope arrays[END_REF]. An de réduire le coût calculatoire, des algorithmes tels que l'EM (Expectation Maximization) ou le SAGE (Space Alternating Generalized Expectation Maximization) ont déjà été proposés [14]. Malheureusement, l'hypothèse de bruit Gaussien est souvent considéré et n'est généralement pas vérié en pratique. Cette non-gaussianité du bruit s'explique par la présence de sources inconnues de très faible intensité ou des interférences causées par des évènements ponctuels ou encore des brouillages articiels, créant ainsi aléatoirement des valeurs aberrantes dans les données mesurées [START_REF] Raza | Spatial ltering of RF interference in radio astronomy[END_REF]. À ce jour et à notre connaissance, la seule tentative de calibration robuste repose sur une modélisation bien spécique du bruit : la loi de Student, avec des variables xxiii indépendantes et identiquement distribuées [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF]. Nous proposons ici d'avoir recours à une plus large gamme de distributions, qui inclut entre autres la loi de Student et le cas classique Gaussien.

Plus particulièrement, nous nous intéressons au processus dit Gaussien composé, qui s'écrit comme le produit d'une variable aléatoire positive, appelée texture, et d'un processus aléatoire Gaussien, le speckle [START_REF] Ollila | Complex elliptically symmetric distributions survey, new results and applications[END_REF][START_REF] Jay | Détection en environnement non gaussien[END_REF]. À partir de ce modèle, nous développons la méthode du maximum de vraisemblance pour estimer à la fois les paramètres de nuisance (texture et speckle) ainsi que les paramètres d'intérêt (matrices de Jones, introduites ci-après). Dans un premier temps, nous nous focalisons sur le cas mono-fréquence.

Considérons D sources de calibration et M antennes dans le réseau de capteurs. Chaque signal incident se décompose selon deux axes de polarisation orthogonales (x, y). Le signal émis par la i-ème source e i et la tension mesurée par la p-ème antenne v i,p (θ) sont liés selon la relation suivante [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF][START_REF] Ker | Radio AGN evolution with low frequency radio surveys[END_REF] v i,p (θ) = J i,p (θ)e i (1) où J i,p (θ) est une matrice de Jones de taille 2 × 2, paramétrée par le vecteur inconnu θ. Elle tient compte de toutes les perturbations introduites le long du trajet depuis la source i jusqu'à l'antenne p. Puisqu'il existe une matrice de Jones diérente pour chaque chemin i-p, le nombre total de matrices à estimer est donc DM .

Le principe d'un interféromètre est de mesurer les corrélations des signaux récoltés par chaque paire d'antennes (p, q). Dans le cas non bruité, cette corrélation s'écrit comme suit

V pq (θ) = E D i=1 v i,p (θ)v H i,q (θ) = D i=1 J i,p (θ)C i J H i,q (θ) avec p < q, p, q ∈ {1, . . . , M } 2 , (2)
où C i = E{e i e H i } décrit l'état de polarisation de la i-ème source, supposé connu, si il s'agit d'une source de calibration. De manière équivalente, nous pouvons considérer

ṽpq (θ) = vec V pq (θ) = D i=1 s i,pq (θ) (3)
où s i,pq (θ) = J * i,q (θ) ⊗ J i,p (θ) c i et c i = vec(C i ).
En pratique, les données sont perturbées par un vecteur bruit n pq qui prend en compte le bruit ambiant Gaussien mais aussi la présence d'éventuelles valeurs aberrantes. Ces dernières peuvent être dues à des phénomènes d'interférences ou encore à la présence de sources de très faible intensité qui nous sont inconnues et inaccessibles dans un premier temps.

Finalement, le vecteur de corrélation bruité de taille 4 × 1 s'écrit

v pq = ṽpq (θ) + n pq . (4) En notant B = M (M -1)
2 le nombre total de paires d'antennes (p, q), le vecteur global de données y ∈ C 4B×1 , regroupant toutes les corrélations mesurées par l'interféromètre, est donné par

y = [v T 12 , v T 13 , . . . , v T (M -1)M ] T = D i=1 s i (θ) + n (5)
tel que s i (θ) = s T i,12 (θ), s T i,13 (θ), . . . , s T i,(M -1)M (θ)

n pq = √ τ pq g pq , (6)
où le facteur de puissance réel τ pq est positif et aléatoire alors que le vecteur g pq suit une loi gaussienne de moyenne nulle, i.e., g pq ∼ CN (0, Ω). Pour éviter tout problème d'ambiguïté, une contrainte, purement arbitraire, est requise sur la matrice de covariance Ω. Dans ce travail, nous supposerons que tr {Ω} = 1.

En supposant des mesures indépendantes entre elles, la vraisemblance admet l'expression suivante

p(y|θ, τ , Ω) = pq 1 |πτ pq Ω| exp - 1 τ pq u H pq (θ)Ω -1 u pq (θ) , (7)
où le vecteur τ = [τ 12 , τ 13 , . . . , τ (M -1)M ] T regroupe les réalisations de la texture et u pq (θ) = v pq -ṽpq (θ). Ainsi, la log-vraisemblance s'écrit

ln p(y|θ, τ , Ω) = -4B ln π -4 pq ln τ pq -B ln |Ω| - pq 1 τ pq u H pq (θ)Ω -1 u pq (θ). (8)
Le principe de la méthode proposée consiste à estimer séquentiellement les paramètres θ, τ et Ω. Ces estimations sont obtenues en maximisant la log-vraisemblance (8) par rapport à chaque paramètre inconnu, les autres paramètres étant supposés xés [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF]. Nous remarquons qu'il est nécessaire de spécier la densité de probabilité du paramètre texture τ pq dans [START_REF] Jongerius | An end-to-end computing model for the Square Kilometre Array[END_REF]. Pour cela, nous pouvons considérer diérentes distributions a priori telles que les lois Gamma, inverse-Gamma, exponentielle et inverse-gaussienne, qui fournissent des formes spéciques pour le modèle Gaussien composé et permettent de dériver diérents estimateurs basés sur le MAP (Maximum A Posteriori) [START_REF] Van Der Tol | Bayesian estimation for ionospheric calibration in radio astronomy[END_REF]. Cependant, supposer des paramètres τ pq inconnus et déterministes permet d'assurer plus de exibilité et de robustesse [START_REF] Conte | Recursive estimation of the covariance matrix of a compound-gaussian process and its application to adaptive CFAR detection[END_REF][START_REF] Pascal | Covariance structure maximum-likelihood estimates in compound gaussian noise: Existence and algorithm analysis[END_REF]. Avec cette hypothèse et la méthodologie itérative adoptée, nous obtenons les estimations des paramètres inconnus décrites ci-dessous 1) Estimation de τ pq : Annuler la dérivée de [START_REF] Dewdney | The Square Kilometre Array[END_REF] par rapport à τ pq conduit à

τ pq = 1 4 u H pq (θ)Ω -1 u pq (θ). (9)
2) Estimation de Ω : Par une méthodologie similaire, nous obtenons l'expression suivante pour la partie speckle

Ωt+1 = 4 B pq u pq (θ)u H pq (θ) u H pq (θ) Ωt -1 u pq (θ) (10

Matrices de Jones non structurées

Pour réaliser la minimisation en [START_REF]Multisource self-calibration for sensor arrays[END_REF], des méthodes avec un coût de calcul réduit sont requises.

En eet, le vecteur d'intérêt θ contient un très grand nombre d'éléments. Notons que ce dernier peut se décomposer comme suit

θ = [θ T 1 , . . . , θ T D ] T = [θ T 1,1 , . . . , θ T 1,M , . . . , θ T D,1 , . . . , θ T D,M ] T , (13)
où θ i,p ∈ R 8×1 décrit le trajet i-p. Ainsi, il est judicieux d'écrire J i,p (θ) ⇒ J i,p (θ i,p ).
Pour réduire les temps de calcul, nous avons recours à l'algorithme EM qui calcule les estimations des paramètres d'intérêt à l'aide de deux étapes successives : les étapes E et M, et tend asymptotiquement vers l'estimation au sens du maximum de vraisemblance [START_REF] Moon | The Expectation-Maximization algorithm[END_REF][START_REF] Feder | Parameter estimation of superimposed signals using the EM algorithm[END_REF]. Dans notre cas, l'étape E se résume à calculer l'espérance conditionnelle de données dites complètes, sachant les données observées et les paramètres estimés. Ensuite, l'étape M maximise la vraisemblance conditionnelle, soit de manière numérique avec l'algorithme LM (Levenberg-Marquardt ) par exemple [START_REF] Nocedal | Numerical optimization[END_REF], soit de manière analytique si possible. Ainsi, [START_REF]Multisource self-calibration for sensor arrays[END_REF] peut se décomposer en plusieurs sous-problèmes de plus faible complexité : l'optimisation se fait par rapport à θ i ∈ C 4M ×1 au lieu de θ ∈ C 4DM ×1 . Les deux étapes successives de l'EM s'écrivent comme suit 1) Étape E : Soit le vecteur de données complètes w = [w T 1 , . . . , w T D ] T , où

w i = s i (θ i ) + n i (14) tel que y = D i=1 w i et n = D i=1 n i . Les vecteurs bruit n i sont supposés statistiquement indépendants et générés selon n i ∼ CN (0, β i Ψ) avec D i=1 β i = 1 et Ψ = bdiag τ 12 Ω, . . . , τ (M -1)M Ω . (15)
L'espérance conditionnelle caractéristique de l'étape E de l'EM est alors donnée par

ŵi = E w i |y; θ, τ , Ω = s i (θ i ) + β i y - D l=1 s l (θ l ) . (16)
2) Étape M : Une fois les w i obtenus avec (16) pour i ∈ {1, . . . , D}, les θ i sont estimés en maximisant la vraisemblance suivante

p( ŵ|θ, τ , Ω) = D i=1 1 |πβ i Ψ| exp -ŵi -s i (θ i ) H (β i Ψ) -1 ŵi -s i (θ i ) . ( 17
)
xxvi Il faut donc considérer la fonction de coût φ i (θ

i ) = ŵi -s i (θ i ) H (β i Ψ) -1 ŵi -s i (θ i ) pour la i-ème source.
Auparavant, nous avons chercher à réaliser la minimisation en [START_REF]Multisource self-calibration for sensor arrays[END_REF] pour chaque θ i . De manière similaire, nous allons par la suite maximiser la vraisemblance en (17) pour chaque θ i,p , de manière alternée, grâce à l'algorithme BCD (Block Coordinate Descent) [START_REF] Friedman | Pathwise coordinate optimization[END_REF] et ce, toujours dans un souci de réduction de complexité.

Nous cherchons à minimiser la fonction φ i (•) par rapport à θ i,p ∈ C 4×1 , tout en xant les autres θ i,q pour q = p. La procédure est répétée séquentiellement pour p ∈ {1, . . . , M } et fournit

les résultats suivants θi,p =    (Σ H i A i,p Σ i + Υ H i Ãi,p Υ i ) -1 (Σ H i A i,p w i,p + Υ H i Ãi,p wi,p ) pour 1 < p < M (Σ H i A i,p Σ i ) -1 Σ H i A i,p w i,p pour p = 1 (Υ H i Ãi,p Υ i ) -1 Υ H i Ãi,p wi,p pour p = M (18) où • w i,p = [w T i,p(p+1) , . . . , w T i,pM ] T et wi,p = [w * T i,1p , . . . , w * T i,(p-1)p ] T , • Σ i = [Σ T i,p+1 , • • • , Σ T i,M ] T et Υ i = [Υ * T i,1 , • • • , Υ * T i,p-1 ] T , • A i,p = bdiag{β i τ p(p+1) Ω, . . . , β i τ pM Ω} -1 et Ãi,p = bdiag{β i τ 1p Ω * , . . . , β i τ (p-1)p Ω * } -1 .
Nous avons également utilisé les notations suivantes :

Σ i,q =     α i,q β i,q 0 0 0 0 α i,q β i,q γ i,q ρ i,q 0 0 0 0 γ i,q ρ i,q     , Υ i,q =     λ i,q µ i,q 0 0 ν i,q ξ i,q 0 0 0 0 λ i,q µ i,q 0 0 ν i,q ξ i,q     avec • α i,q = q * i 1 c i 1 + q * i 2 c i 3 , β i,q = q * i 1 c i 2 + q * i 2 c i 4 , γ i,q = q * i 3 c i 1 + q * i 4 c i 3 , ρ i,q = q * i 3 c i 2 + q * i 4 c i 4 , • λ i,q = q i 1 c i 1 + q i 2 c i 2 , µ i,q = q i 1 c i 3 + q i 2 c i 4 , ν i,q = q i 3 c i 1 + q i 4 c i 2 , ξ i,q = q i 3 c i 3 + q i 4 c i 4 , • c i = [c i 1 , c i 2 , c i 3 , c i 4 ] T , J i,p (θ i,p ) = p i 1 p i 2 p i 3 p i 4 , J i,q (θ i,q ) = q i 1 q i 2 q i 3 q i 4 , • θ i,p = [p i 1 , p i 2 , p i 3 , p i 4 ] T et θ i,q = [q i 1 , q i 2 , q i 3 , q i 4 ] T .
L'algorithme de calibration robuste que nous proposons est résumé ci-dessous.

Algorithme: Calibration robuste basée sur un modèle Gaussien composé entrée : D, M , B, C i , β i , y sortie :

θ initialisation: Ω ← Ω init , τ ← τ init , θ ← θ init while critère d'

Matrices de Jones structurées

Une matrice de Jones peut se décomposer comme le produit de plusieurs matrices, chacune représentant un eet physique bien particulier. Dans le cas d'un régime de calibration appelé régime 3 [START_REF] Lonsdale | Calibration approaches, LFD memo 015[END_REF], exposé en Figure 2, les antennes du réseau sont relativement proches les unes des autres et admettent chacune un très large champ de vision, comme c'est le cas à l'échelle d'une station du LOFAR ou du SKA [START_REF] Wijnholds | Calibration challenges for future radio telescopes[END_REF]. Dans un tel contexte, nous considérons la décomposition suivante pour les matrices de Jones [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF][START_REF] Thompson | Interferometry and synthesis in radio astronomy[END_REF][START_REF] Yatawatta | Reduced ambiguity calibration for LOFAR[END_REF]]

J i,p (θ i,p ) = G p (g p )H i,p Z i (ϕ i )F i (ϑ i ) ( 19
) où i ∈ {1, . . . , D}, p ∈ {1, . . . , M } et θ i,p = [ϑ i , ϕ i , g T p ]
T . La seule matrice supposée connue ici est H i,p , qui tient compte du diagramme de rayonnement des antennes ainsi que du retard de propagation géométrique. Elle peut être supposée connue grâce à des simulations électromagnétiques eectuées au préalable et à l'information disponible fournie par la position de l'antenne et la direction de visée [START_REF] Hamaker | Understanding radio polarimetry[END_REF][START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF][START_REF] Kazemi | Blind calibration for radio interferometry using convex optimization[END_REF]. Du fait de la structure compacte du réseau d'antennes considéré, la même partie de l'ionosphère est observée et les eets environnementaux sont supposés identiques pour l'ensemble des antennes. Plus précisement, le passage dans l'ionosphère entraîne un déphasage qui s'écrit ainsi qu'une rotation du plan de polarisation, appelée rotation de Faraday et donnée par

Z i (ϕ i ) = exp jϕ i I 2
F i (ϑ i ) = cos(ϑ i ) -sin(ϑ i ) sin(ϑ i ) cos(ϑ i ) . (21)
Enn, chaque antenne est décrite par un gain instrumental caractéristique de la chaîne de réception, tel que

G p (g p ) = diag{g p }. (22)
Sur la Figure 2, nous exposons le principe du régime 3 : la même portion de l'ionosphère est considérée mais une multitude de sources peuvent être visibles dans le champ de vision. Notons que dans le cas où les antennes sont relativement éloignées les unes des autres et admettent toujours chacune un très large champ de vision, comme c'est le cas à l'échelle des plusieurs stations, il s'agit alors du régime 4 où nous pouvons introduire un déphasage ionosphérique ϕ i,p et un angle de rotation ϑ i,p pour chaque source i ∈ {1, . . . , D} et chaque antenne p ∈ {1, . . . , M }.

Après avoir estimé les entrées des matrices de Jones à l'aide de l'algorithme décrit dans le cas non structuré, il est possible de déduire les paramètres physiques qui interviennent dans [START_REF] Ker | Radio AGN evolution with low frequency radio surveys[END_REF].

Pour cela, nous considérons une fonction coût de type moindres carrés et optimisons par rapport à chacun des paramètres inconnus de manière alternée. Ainsi, nous obtenons les estimations suivantes 1) Estimation de g p : Après calculs, chaque gain est donné par

[ĝ p ] k = D i=1 [W * i,p ] k,k -1 D i=1 [X * i,p ] k,k (23)
xxix

où X i,p = R i,p ĴH i,p , W i,p = R i,p R H i,p et R i,p = H i,p Z i F i .
2) Estimation de ϕ i : Après dérivation et égalisation à zéro, nous obtenons le résultat suivant exp 2j φi =

Tr M i,p Tr M H i,p (24)
où

M i,p = M p=1 Ĵi,p F H i H H i,p G H p et φi peut directement en être déduit.
3) Estimation de ϑ i : Une procédure numérique pour une dimension est nécessaire pour l'estimation des angles de Faraday. Cette dernière peut facilement être réalisée en un temps limité à l'aide d'algorithmes classiques de type Newton ou descente de gradient.

Pour le moment, nous nous sommes restreints au cas mono-fréquence. Or, en pratique, les systèmes d'observation actuels exploitent une large gamme de fréquences. Il nous semble donc essentiel de prendre en compte l'aspect multi-fréquentiel.

Méthode de calibration robuste multi-fréquentielle

Pour illustrer l'utilisation de notre algorithme dans un contexte multi-fréquentiel, nous choisissons ici de considérer le régime 4, où les stations d'antennes sont éloignées les unes des autres, et la décomposition en [START_REF] Ker | Radio AGN evolution with low frequency radio surveys[END_REF] peut s'écrire [25, 33]. Notons que l'adaptation de l'algorithme multi-fréquentiel au régime 3 peut directement être déduite à partir de l'étude proposée. Pour les gains électroniques g [f ] , nous pouvons supposer une certaine régularité dans la variation fréquentielle à l'aide d'un modèle polynomial. Ainsi, nous introduisons

J [f ] i,p (θ [f ] ) = G [f ] p (g [f ] p )H [f ] i,p Z [f ] i,p (ϕ [f ] i,p )F [f ] i,p (ϑ [f ] i,p ) (25) où ϕ
[f ] i,p ∝ 1 f et ϑ [f ] i,p ∝ 1 f 2
g [f ] p = B[f] zp (26) où B[f] = b[f] T ⊗ I 2 et [ b[f] ] k = f -f 0 f 0 k-1 (27)
avec k ∈ {1, . . . , Ñ } et une fréquence de référence f 0 . L'ordre Ñ est choisi de manière arbitraire et les variables latentes zp nécessitent d'être estimées.

En supposant des mesures indépendantes entre paires d'antennes mais aussi entre fréquences et le terme de pénalisation comme suit

f ∈ F = {f 1 , . . . , f F },
L {θ [f ] } f ∈F , z, {x [f ] } f ∈F = f ∈F L [f ] θ [f ] , z, x [f ] (28) tel que L [f ] θ [f ] , z, x [f ] = l [f ] θ [f ] + h [f ] θ [f ] , z, x [f ] , où le terme d'attache aux données s'écrit l [f ] (θ [f ] ) = pq 1 τ [f ] pq u [f ] H pq (θ [f ] )Ω [f ] -1 u [f ] pq (θ [f ] ) (29)
h [f ] θ [f ] , z, x [f ] = 2 x [f ] H θ [f ] -B [f ] z + ρ||θ [f ] -B [f ] z|| 2 2 . ( 30
)
Le vecteur x [f ] fait référence aux variables duales du problème pour une fréquence f , z inclut l'ensemble des variables cachées/latentes tandis que la matrice de couplage B [f ] tient compte des variations supposées connues des paramètres d'intérêt selon la fréquence. Pour résoudre un tel problème séparable en fréquence, nous utilisons une méthode dite distribuée qui permet de répartir la charge de calcul massive sur plusieurs agents [START_REF] Bertsekas | Parallel and distributed computation: numerical methods[END_REF][START_REF] Yatawatta | Distributed radio interferometric calibration[END_REF]. Ainsi, la calibration est réalisée localement au-niveau de chaque agent pour chaque fréquence tandis qu'un centre de fusion de l'information permet d'assurer la communication entre les diérents agents, d'où une meilleure calibration possible puisque l'ensemble des fréquences est pris en compte et non pas une fréquence après l'autre de manière indépendante. Pour illustrer cela, nous montrons sur la Figure 3 l'évolution de l'EQM en fonction du RSB pour un paramètre représentatif des matrices de Jones dans deux cas spéciques : i) calibration multi-fréquentielle avec une procédure distribuée et ii) calibration mono-fréquentielle où chaque fréquence est traitée de manière indépendante.

Une méthode appropriée pour résoudre un problème d'optimisation sous contraintes de manière distribuée est l'algorithme ADMM (Alternating Direction Method of Multipliers) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] qui se compose de trois étapes successives, répétées itérativement,

• θ[f] t+1 = argmin θ [f ] L [f ] θ [f ] , (ẑ) t , x[f] t réalisée localement par chaque agent (31)
• (ẑ) t+1 = argmin

z f ∈F h [f ] θ[f] t+1 , z, x[f] t réalisée globalement au centre de regroupement (32) • x[f] t+1 = x[f] t +ρ θ[f] t+1 -B [f ]
(ẑ) t+1 réalisée localement par chaque agent [START_REF] Van Der Tol | Ionospheric calibration for the LOFAR radio telescope[END_REF] Notons que l'étape [START_REF] Kazemi | Blind calibration for radio interferometry using convex optimization[END_REF] donne directement une expression analytique tandis que l'étape [START_REF] Yatawatta | Reduced ambiguity calibration for LOFAR[END_REF] nécessite une approche itérative en minimisant de manière alternée le Lagrangien L à chaque ϑ

[f ] par rapport xxxi (a) (b) (c)
[f ] i,p , ϕ [f ] i,p et g [f ]
p . Après calculs, des expressions analytiques peuvent être obtenues pour les gains électroniques ainsi que les déphasages ionosphériques. En ce qui concerne les angles de Faraday, une optimisation numérique est requise. Par souci de clarté, les expressions analytiques ainsi que le détail de l'optimisation alternée pour résoudre (31) sont donnés en section 4.3.2.

Simulations réalistes sous MeqTrees

Ce travail a aussi été l'occasion d'utiliser l'outil de simulation MeqTrees qui permet de générer des données interférométriques telles que celles mesurées par de réels instruments astronomiques [START_REF] Noordam | The MeqTrees software system and its use for thirdgeneration calibration of radio interferometers[END_REF] For sake of clarity and coherence, the work presented in this document does not account for the study made for separated arrays and structured noise covariance matrix [C1] as well as the derivation of Bayesian lower bounds in the Random Matrix Theory (RMT) framework [C3]. In order to reach the theoretical optimal performances of this new generation of radio interferometers, signal processing challenges must be overcome in terms of data correction and image synthesis. Indeed, estimation of all perturbation eects along the radio signal propagation path, namely the calibration process, is of critical importance in order to reconstruct images with no distortions. In this manuscript, we focus on the develoment of signal processing calibration algorithms, as it is a cornerstone of the imaging step in radio astronomy. To the best of our knowledge, the only alternative robust algorithm using a non-Gaussian noise assumption is presented in [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF][START_REF] Kazemi | Robust radio interferometric calibration using the tdistribution[END_REF] and based on the Student's t noise modeling. Nevertheless, this modeling may not be the optimal way to incorporate the contribution of outliers in radio data as the distribution of such data is not known in practical scenario.

• Propose estimation algorithms with low computational complexity and ecient implementation.

• Study two major scenarios: calibration of interferometric systems with wide eld-of-view and closely located (one station) or largely separated (multiple stations) radio antennas.

These two dierent scales of observation refer, respectively, to calibration regimes 3 and 4, exposed in section 2.2.3 and are the most adapted for recent and next generation radio interferometers [START_REF] Lonsdale | Calibration approaches, LFD memo 015[END_REF].

• Exploit the multi-frequency case where a wide frequency range is considered and ecient processing of data among all dierent frequencies is required.

• Characterize the estimation performances thanks to statistical studies and realistic reconstructed images in which weak unknown background sources can be revealed.

Achieved results

In this document, we propose robust calibration algorithms which are iterative, based on the Maximum Likelihood (ML) estimator and the Compound-Gaussian (CG) noise modeling. The CG distribution encompasses a wide range of distributions including the typical Gaussian, the Student's t, the Laplace and the Cauchy cases. Therefore, we choose to consider such exible noise modeling to achieve robustness. In order to reduce the computational cost, eorts are being made to obtain closed-form expressions for all estimates when possible, by combining the use of the Expectation Maximization (EM) and Block Coordinate Descent (BCD) algorithms.

In radio astronomy, calibration amounts to estimate specic 2 × 2 matrices, called Jones matrices as described in section 2.3, which stand for all the introduced environmental and instrumental perturbation eects [START_REF] Hamaker | Understanding radio polarimetry[END_REF][START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF][START_REF] Thompson | Interferometry and synthesis in radio astronomy[END_REF]. The rst proposed robust algorithm in this document, namely the Relaxed Concentrated ML Estimator (RCMLE), operates in a mono-frequency scenario and is adapted for both non-structured and structured cases. This refers, respectively, to estimation of all entries of all Jones matrices and estimation of specic parameters in physical models (regime 3 as a rst step). Statistical analysis with the Mean Square Error (MSE) and the Cramer-Rao Bound (CRB) highlight that the proposed scheme leads to more accurate estimation, thus, more robustness w.r.t the presence of weak unknown sources in the background and interferences. Let us note that the designed RCMLE is based on a relaxed version of the ML, i.e., a deterministic assumption but its alternative, namely the Bayesian approach, is proposed as well in this document.

MANUSCRIPT STRUCTURE 3

To extend our study to the multi-frequency case, we adopt a distributed approach with multiple local agents and one fusion center, unlike the classical strategy which uses a single centralized processor. We exploit the specic varition of parameters across frequency and conduct the study for the structured case. This leads to the proposed Multi-frequency Robust Calibration Algorithm (MRCA) designed for regime 4 with a possible straightforward adaptation to regime 3. Numerical simulations compare the estimation performances w.r.t. state-of-the art non-robust and/or mono-frequency cases (per-channel calibration).

Finally, to complete the analysis, we make use of the software MeqTrees for simulation of realistic radio data and visual indication of our improved performances in terms of calibration and image reconstruction [START_REF] Noordam | The MeqTrees software system and its use for thirdgeneration calibration of radio interferometers[END_REF].

Manuscript structure

The presented document relies mainly on the four following chapters

• Chapter 2 is dedicated to state-of-the-art calibration algorithms. We rst present the context under consideration with a brief operating principle of interferometry and the new challenges at stake accompanying the development of advanced radio interferometers. The data model supported by Jones matrices is also introduced and most well-known calibration strategies are highlighted. A short description of the imaging procedure completes the chapter.

• Chapter 3 develops the proposed robust calibration algorithm based on th CG noise modeling. As a rst step, the frequency dependence is not considered in this chapter.

• Chapter 4 adapts our methodological tools to the multi-frequency case. A distributed strategy is employed and robustness is still ensured.

• Lastly, chapter 5 introduces realistic data simulation with the software MeqTrees which helps us to visualize reconstructed residual images and reveal the presence of weak noncalibration sources, hidden in the noisy background.

Publications

The work described in this document has led to the following publications Peer-reviewed international journals With a single telescope to probe the sky, the corresponding angular resolution is given by

θ = λ D (2.1)
where λ is the incident wavelength and D is the diameter of the instrument. Thus, a better resolution requires a larger diameter, which can involve cost and construction complexity issues.

To overcome these problems, the concept of interferometry was applied by considering an array with multiple antennas instead of one. The SKA [START_REF] Jongerius | An end-to-end computing model for the Square Kilometre Array[END_REF][START_REF] Dewdney | The Square Kilometre Array[END_REF][START_REF] Ferrari | French SKA White Book -The French community towards the Square Kilometre Array[END_REF] is the upcoming new generation of radio antennas, which will operate in a large observation bandwidth (up to tens of GHz) with a total collecting area of roughly one square kilometre. Expected to be the most sensitive and fastest astronomical instrument to probe the radio sky, its stations will be distributed in South Africa and Australia, composed of both traditional dishes and small dipoles. The construction is set to begin in 2018 and will be extended to its full size until 2030, in two consecutive phases. Precursor facilities of the SKA are already currently running: the Australian Square Kilometre Array Pathnder (ASKAP), the Murchison Wideeld Array (MWA) in Australia and the Karoo Array Telescope (MeerKAT) in South Africa.

Principles of interferometry

As mentioned before, interferometry consists in measuring correlations, called visibilities in radio astronomy, for two specic antennas, specically distributed over an area at a given wavelength [START_REF] Thompson | Interferometry and synthesis in radio astronomy[END_REF]. Due to Earth's rotation, baseline distances and orientations w.r.t. the sky eld vary over time leading to more data points and high-resolution synthesis imaging. To give an idea, the longest baselines of the SKA will reach 3000 km at least, while the LOFAR attains 1000 km and the VLA, 36 km. By denition, a correlation V is given by the following expectation [START_REF] Ker | Radio AGN evolution with low frequency radio surveys[END_REF] V (r 1 , r 2 ) = E e(r 1 )e * (r 2 )

(2.3)
where, for sake of clarity, the electrif eld e is assumed scalar and unpolarized, measured for two points in space, r 1 and r 2 . For two locations, separated by a given propagation delay τ , one obtains [START_REF] Thompson | Interferometry and synthesis in radio astronomy[END_REF] V

(r 1 , r 2 , τ ) = lim T →∞ 1 2T T -T e(r 1 , t)e * (r 2 , t -τ )dt.
(2.4)

In a suitable coordinate system, it is possible to link this spatial correlation with the source intensity. The relative distance r 1 -r 2 is usually expressed in units of wavelength, given by three coordinates (u, v, w), in which w is omitted if the coplanar approximation holds or the eld-ofview is narrow. In the two-dimensional case, the visibility is given by the Fourier Transform of the source intensity, i.e., [START_REF] Levanda | Synthetic aperture radio telescopes[END_REF] V (u, v) = I(l, m)e -2πj(ul+vm) dldm

(2.5)

CONTEXT OF STUDY 9

where (l, m) refers to the angular coordinates of the source (right ascension (RA) and declination (DEC) for exemple), I(l, m) is the intensity distribution at a given point and V (u, v) is the measured correlation for a given baseline. In the ideal case, when the entire (u, v) plane is sampled by visibility measurements, the source intensity is given by the following Fourier inversion ul+vm) dudv.

I(l, m) = V (u, v)e 2πj(
(2.6)

The mathematical relationship in (2.5) is a particular expression of the fundamental Van Cittert-Zernike theorem [26, chp. 14], obtained by Van Cittert and Zernike [START_REF] Zernike | The concept of degree of coherence and its application to optical problems[END_REF], which is based on coherence theory and radiation emission by an incoherent electric eld in the far-eld. This theorem stems from the optical domain [START_REF] Born | Principles of optics -Electromagnetic theory of propagation, interference and diraction of light -seventh[END_REF] and requires specic conditions: the emitting source is spatially incoherent, i.e., the radiation emitted by one point of the source is statistically independent from any other point of the source (which is the case for most of the astronomical sources expect for pulsars and masers). In addition, the angular dimensions of the source must be typically inferior to one degree, the spectral width is quasi monochromatic and the source is in the far-eld zone. It is worth mentioning that this relationship was initially obtained in the case of free-space between source and antennas, implying a uniform refractive index and an undistorted wavefront. Still, the generalization is achievable for an heterogeneous propagation medium and is referred to as Hopkins' formula. The principle of data reconstruction with an interferometric system is exposed in ul+vm) dudv.

I D (l, m) = V (u, v)S(u, v)e 2πj(

Challenges of the new generation of radio interferometers

Advanced phased arrays as the LOFAR and the SKA provide considerable benets such as huge collecting area and large spatial coverage, ne sensitivity and exibility, the possibility to observe in multiple directions simultaneously, a quick response time, ... However, great benets necessarily produce some diculties. Indeed, new challenges are at stake: a huge number of antennas with wide eld-of-view, looking through an heterogeneous propagation medium, leads to signicant direction dependent perturbation eects. Besides, the phased-array observational system has multi-frequency and multi-beam forming capabilities, polarization of waves needs to be considered and the theoretically achievable angular, temporal and spectral resolution is high. Antennas being ominidirectional and particularly sensitive, they can detect an important number of sources in their eld-of-view. And in some cases, the source of interest reveals to be more complex than a simple point source, with constant intensity. Indeed, its structured diuse 2.2. IMPORTANCE OF CALIBRATION 11 emission and its polarization state may change according to the frequency channel and along the signal propagation path.

With wide eld-of-view and well-separated stations, the next generation array cannot be considered coplanar and spherical geometry is required. Indeed, all antennas are not located in a given plane at the same altitude. In this particular case, spatial coherency and brightness distribution are no longer linked thanks to the two-dimensional Fourier Transform in (2.5), since the coordinate w is essential and curvature of the Earth cannot be neglected. To solve these issues, faceting-based methods have been introduced as well as projection methods [START_REF] Mcmullin | CASA architecture and applications[END_REF][START_REF] Myers | CASA synthesis & single dish reduction cookbook[END_REF][START_REF] Cornwell | The noncoplanar baselines eect in radio interferometry: The W-projection algorithm[END_REF].

Last but not least, the new radio instruments are particularly sensitive and the amount of data to process is tremendous. We need to deal with large volume of data depending on the RA, the DEC and the frequency 1 .

Therefore, to meet the theoretical optimal performances, many signal processing challenges need to be overcome, among which calibration, imaging and data processing [START_REF] Wijnholds | Calibration challenges for future radio telescopes[END_REF][START_REF] Noordam | LOFAR calibration challenges[END_REF][START_REF] Wijnholds | Signal processing challenges for radio astronomical arrays[END_REF]. In the following, we focus on the calibration step. Indeed, visibilities have assumed to be ideal so far, i.e., uncorrupted. However, in practice, observed visibilities dier from the true visibilities as instrumentation is not perfect and physical disturbing phenomena are occuring. Thus, correcting these eects is essential in order to produce high dynamic range images with advanced observation systems.

2.2

Importance of calibration

Perturbation eects

When an incoming electromagnetic eld reaches an array, it has been aected by various perturbation eects and physical distortions which involve propagation delay and magnitude loss. Here, we draw a distinction between, on the one hand, instrumental and electronic eects linked to the receiver chain and on the other hand, environmental eects due to propagation within the atmosphere.

• Instrumental eects: each antenna is dened by an unknown electronic gain and anisotropic beam pattern due to direction dependency and selectivity leakage. Composed of a main beam and sidelobes, the beam pattern accounts for sensitivity of antennas to radiation from any viewing direction. Electromagnetic simulations and modeling provide some a priori knowledge about it.

However, because of sidelobes, an antenna is also sensitive in directions dierent from the one of interest. Thus, bright sources can hinder the observation of faint sources, which is particularly troublesome in the case of the LOFAR owing to its relatively large sidelobes. In addition, let us note that the beam pattern diers per station, w.r.t. time and frequency due to Earth's rotation. In the receiver chain, temperature may vary as well and other properties along the signal propagation path are changing through time until reaching the data processing center, thus aecting calibration parameters as gains.

• Environmental eects: the troposphere and the ionosphere are specic layers of the Earth's atmosphere. They entail time delays due to temporally and spatially varying diraction and refraction phenomena, but also eects on the polarization of the signal as Faraday rotation [START_REF] Davies | Ionospheric radio[END_REF] which is a rotation of the polarization plane of the wave while it goes through a charged medium as the atmosphere. This rotation is due to dierent propagation speeds between the two polarized components of the incoming wave. One prevailing eect is thus propagation within the 1 Let us note that astronomical objects are identied by equatorial coordinates: the RA and the DEC, equivalent to celestial coordinates: latitude and longitude.

CHAPTER 2. PROBLEM SETUP AND STATE-OF-THE-ART

ionosphere, the ionized upper part, which leads to signicant rapidly changing phase delays. We specify this point in the following remark.

Remark: The ionosphere is a specic layer of the atmosphere, which introduces time and phase delays, resulting in a shifted position of the source. Thus, the apparent position diers from the exact location. These visible consequences are due to non-uniform refractive index in the propagation medium and spatial, temporal uctuations in the electron density. More specically, the sun emits radiations which induce partial ionization in the atmosphere through injection of charged particles. This is called the solar wind. Therefore, free electrons emerge and slow the progression of electromagnetic waves. The movement of these free electrons is determined by the magnetic eld, ionization, dissociation and ion recombination processes.

We can view the ionosphere as a dynamic, disruptive and heterogeneous environment which introduces phase shifts proportional to the wavelength, thus, inversely proportional to the frequency [START_REF] Van Der Tol | Bayesian estimation for ionospheric calibration in radio astronomy[END_REF][START_REF] Davies | Ionospheric radio[END_REF]. This means that ionospheric perturbations are even more prevailing at low frequencies, especially for the LOFAR and the SKA which operate at low radio frequencies.

Ionospheric disturbances are also temporally varying since at night, the number of electrons drops while it is higher during the day. Typically, at night and at the zenith, the number of electrons is given by 10 16 electrons per square metre, which is equivalent to one unit of Total Electron Content (T EC). The T EC refers to the integral of electron density along the viewing direction and is directly related to the propagation time dierences. Let us mention that correcting for artefacts due to the ionosphere remains a critical step in many scientic elds and not just in radio astronomy [START_REF] Jung | An improvement of ionospheric phase correction by multipleaperture interferometry[END_REF].

Aim of calibration

In radio astronomy, calibration is of the utmost importance before synthesis imaging and it directly deals with visibility measurements which may be aected by the perturbations, listed in section 2.2.1. In Figure 2.3, we show a simulated sky with undisturbed sources of varying intensity on the left hand side [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF]. When the brightest sources suer from perturbations, the resulting image is given on the right hand side. Thus, severe distortions and artefacts appear, making it dicult to reveal faint sources. These distortions are due to the simulated errors along the brightest sources but also to the incomplete Fourier sampling, leading to strong sidelobe patterns which overwhelm the image.

To have a chance of recovering the initial image, unknown perturbation eects need to be estimated and corrected. Afterwards, we can remove the contribution of the bright sources from the data, thanks to tables providing their position and intensity accurately [START_REF] Kimball | Ivezi¢, A unied catalog of radio objects detected by NVSS, FIRST, WENSS, GB6, and SDSS[END_REF][START_REF] Baars | The absolute spectrum of Cas A-an accurate ux density scale and a set of secondary calibrators[END_REF]. In the residual image, only the faintest sources remain visible which is of interest for the astronomers since these sources are the most dicult to access. Let us remark that most sources have an intensity lower than the noise level (except for particularly strong ones like the sun, Cygnus A (CygA) and Cassiopeia A (CasA) radio stars). Getting more knowledge about the sky model, the instrumentation and propagation conditions constitutes the ultimate goal of the calibration process.

Calibration regimes

Calibration of the ionosphere is addressed in dierent ways depending on the type of sensor array. Four dierent scenarios have been dened [START_REF] Lonsdale | Calibration approaches, LFD memo 015[END_REF] and they are presented in Figure 2.4. They are distinguished by the aperture A, the size of the eld-of-view F which is projected on the ionospheric layer and the scale P on which ionospheric perturbations are varying: • In regime 1, the eld-of-view F and the aperture A are very small in comparison with the ionospheric irregularity scale P . Consequently, each antenna sees the same part of the ionosphere so the same phase delay is introduced for each antenna in the array. The eld-of-view being narrow, there is no direction dependency: the propagation conditions are the same for all paths within the eld-of-view.

• In regime 2, the antennas are relatively separated from each other but their eld-of-view remains narrow. When we target a specic direction, propagation conditions dier per antenna.

Nevertheless, the same time delay is introduced for each source within the same eld-of-view.

• Regime 3 refers to compact stations, with relatively closely spaced antennas, as in regime 1. By contrast, the eld-of-view F is wide in comparison with the scale of variation P . Therefore, propagation conditions are direction dependent due to the wide eld-of-view but antennas observe through the same part of the ionosphere, leading to similar perturbations. For instance, this amounts to consider a single station of the LOFAR or a subarray of the MWA.

• The most dicult case is regime 4 where stations are well separated from one another and their eld-of-view are wide. Thus, the ionosphere induces per-antenna direction dependent perturbations.

In brief, ionospheric calibration requires to take into account variable properties, depending on the environment and the telescopes. The most adapted scenarios for the next-generation of radio interferometers are regimes 3 and 4 and dierent advanced approaches using array and signal processing techiques can be considered [5860]. The conventional data model that is studied in radio astronomy is based on the so-called Jones matrices [START_REF] Ker | Radio AGN evolution with low frequency radio surveys[END_REF][START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF]. In the next section, we present this formalism for radio interferometer measurement equation [START_REF] Hamaker | Understanding radio polarimetry[END_REF][START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF].

Data model

Let us consider a signal source with two components, each being associated to one specic polarization direction [START_REF] Salvini | Fast gain calibration in radio astronomy using alternating direction implicit methods: Analysis and applications[END_REF][START_REF]StEFCal an alternating direction implicit method for fast full polarization array calibration[END_REF]. Unlike an unpolarized source which is dened by a single scalar, i.e., its intensity, a polarized source is described by the distribution of the electric eld. Thus, the wave nature of light needs to be considered. A priori information can be given about it thanks to survey observations. To describe the polarization state, we usually introduce the following 4

Stokes parameters [START_REF] Thompson | Interferometry and synthesis in radio astronomy[END_REF]: (I, Q, U, V ) with, respectively, an unpolarized, two linearly polarized and one circularly polarized components. This polarizarion of incident waves depends on position, frequency and can be severely aected by perturbations such as the Faraday rotation [START_REF] Davies | Ionospheric radio[END_REF].

As mentioned before, the signal source is dened along two orthogonal polarization directions (horizontal and vertical or circularly right and left depending on the coordinate system), the third direction being the wave propagation direction. To retrieve both information, each antenna of the array is made up of two receptors (or dipoles), each being sensitive to one specic polarization direction, as shown in Figure 2.5. Thus, all the study is conducted in the plane orthogonal to the wave propagation direction. The incoming radio signal is decomposed as a 2 × 1 vector for the i-th source direction:

e i = [e ix , e iy ] T .
(2.8)

The incident electric eld e i and the measured voltage v i,p (θ) are related in the following way

v i,p (θ) = J i,p (θ)e i (2.9)
where the 2 × 2 Jones matrix, noted J i,p (θ), stands for all physical perturbations along the signal propagation path i-p, i.e., from the i-th emitting source to the p-th receiving antenna, Figure 2.5: Interferometric system (after [START_REF] Hamaker | Understanding radio polarimetry[END_REF]).

and is parametrized by unknown vector θ. The relationship is given by a multiplication as propagation of the wave, and thus all undergone transformations, are assumed linear. We note

v i,p (θ) = [v i,px (θ), v i,py (θ)]
T the generated voltage at the p-th dual polarized antenna and we can associate one particular Jones matrix J i,p with each signal propagation path i-p.

The principle of an interferometer is to measure correlations between two specic antennas p and q. The resulting 2 × 2 ideal correlation matrix is given by

V i,pq (θ) = E v i,p (θ)v H i,q (θ) = J i,p (θ)C i J H i,q (θ) (2.10)
where

C i = E e i e H i = I i + Q i U i + jV i U i -jV i I i -Q i
stands for the intrinsic coherency or brightness of the i-th source, describing its polarization state. It is a function of previously introduced Stokes parameters, themselves functions of complex amplitudes of the electric eld in a specic basis [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF]. To take Jones matrices out of the expectation value in (2.10), they are assumed constant during the integration time interval.

Let us note that an alternative formulation of (2.10) exists, which makes use of the 4 × 4

Mueller matrices [50, p. 417], [START_REF] Rau | Advances in calibration and imaging techniques in radio interferometry[END_REF]. Instead of considering the matrix product between two voltage vectors, the outer product is used, leading to the following correlation vector

v i,pq (θ) = J i,p (θ) ⊗ J H i,q (θ) e i ⊗ e H i = J i,p (θ) ⊗ J H i,q (θ)     I i + Q i U i + jV i U i -jV i I i -Q i     (2.11)
in which J i,p (θ) ⊗ J H i,q (θ) describes the combine perturbation eects along the i-th signal path to antennas p and q. Such matrices are notably considered in [START_REF] Tasse | Applying full polarization A-Projection to very wide eld of view instruments: An imager for LOFAR[END_REF] where the A-Projection deconvolution algorithm [START_REF] Bhatnagar | Wide-eld wide-band interferometric imaging: The WB A-projection and hybrid algorithms[END_REF] is able to reconstruct an image, taking into account directional eects CHAPTER 2. PROBLEM SETUP AND STATE-OF-THE-ART as the frequency, time and polarization dependent beam pattern of antennas [START_REF] Bhatnagar | Calibration and imaging challenges at low radio frequencies: An overview of the state of the art, The Low Frequency Radio Universe[END_REF], and correcting for adverse eects. This particular 4 × 4 form is mathematically equivalent to the 2 × 2 one but is more adapted to depict imaging problems as it emphasizes the linear operations performed by the interferometer on the sky distribution [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF]. In our study, we focus on the calibration procedure and the 2 × 2 formalism is thus favoured as it provides a more transparent, more intuitive and less complicated description of the underlying physics behind the signal propagation and it is also more suitable due to computational issues.

When the sky is composed of D multiple sources, assumed uncorrelated and far away, the emitted waves are considered planar when reaching the Earth. Thus, contributions of all sources are linearly added as follows

V pq = D i=1 J i,p (θ)C i J H i,q (θ).
(2.12)

In a more realistic scenario, the contribution of calibration sources in the sky is not written as a discrete sum but rather like a diuse radiation varying in direction. However, in this calibration study, we assume point sources as a rst step. Besides, in practical scenario, we usually add an uncorrelated Gaussian noise N pq in (2.12).

A given Jones matrix J i,p , for one specic propagation path i-p, can be decomposed into the multiplication of many individual terms, each being associated with one particular eect [START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF]:

J i,p = J i,p,1 J i,p,2 • • • J i,p,T . (2.13)
This is what we call a Jones chain [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF], with T assumed dierent eects here. The arrangement is particularly important since it corresponds to the physical order in which eects appear along the path. Thus, the leftmost matrix J i,p,1 stands for the perturbation introduced belatedly, i.e., close to the receiver, as the electronic gains, while the matrices on the right describe the rst introduced perturbations close to the source emission, as the Faraday rotation for instance.

Various physical and geometrical eects can be described by Jones matrices as ionospheric phase shifts due to refraction phenomena but also rotation of the so-called parallactic angle between the reception system and the incident eld or even change of coordinate system which requires transformation matrices (conversion into circular coordinates for instance can reveal to be more adapted in some cases and this choice depends on the specic design of telescopes). Let us also mention the possible interferences between antennas (since components of the polarized electromagnetic signal are not always well separated and one receptor may collect part of the radiation intended for the other receptor), the complex electronic gains of antennas, the lter bandwidth in the electronic chain and the errors due to the correlator [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF][START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF].

In the literature, Jones matrices have their own specic structures [START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF][START_REF] Sault | The Hamaker-Bregman-Sault measurement equation[END_REF]. Let us start with a general decomposition [START_REF] Ker | Radio AGN evolution with low frequency radio surveys[END_REF] J

i,p = G p B p D p E i,p P i,p K i,p T i,p Z i,p F i,p (2.14)
in which G p refers to the per-antenna diagonal electronic complex gain matrix and B p is the electronic frequency bandpass. The direction independent matrix D p stands for cross-leakage and interferences between receptors while direction dependent perturbations include the primary beam E i,p , a projection matrix P i,p which takes into account the orientation of receptors w.r.t.

the sky, the typical geometric pathlength dierence K i,p , the tropospheric phase delay T i,p , the ionospheric phase delay Z i,p and the ionospheric Faraday rotation F i,p . Particular eects represented by G p and E i,p are studied in [START_REF] Yatawatta | LOFAR beamshapes and their use in calibration and imaging[END_REF], where geometry of dipoles is taken into account within the LOFAR and a possible model of the beam pattern is needed. Besides, E i,p , Z i,p and F i,p are considered in [START_REF] Kazemi | Clustered calibration: an improvement to radio interferometric direction-dependent self-calibration[END_REF]. In section 3.2.2, we will introduce the matrix H i,p to account for all assumed known eects, including matrices E i,p and K i,p . In the following, a detailed description of some of the above matrices with simple representations is given

• The matrix K i,p expresses the propagation delay, i.e., a phase shift due to geometrical path dierence. Such phase shifts can be compensated by introducing additional electronic delays in the corresponding directions. This matrix which is proportional to identity, namely a scalar matrix, aects each component of the signal in the same way so there is independency w.r.t.

polarization. Its expression is given by

K i,p = exp -2πj(u p l i + v p m i + w p (n i -1)) I 2 (2.15)
where the position of the p-th antenna is expressed in units of wavelength (u p , v p , w p ) and the direction of the i-th observed source is dened by three coordinates (l i , m i , n i ). Scalar matrices have the same representation in any coordinate system and can be switched to any place of the Jones chain.

• Another typical scalar Jones matrix expresses the phase shift ϕ i,p induced by the ionosphere.

It is written as

Z i,p = exp jϕ i,p I 2 .
(2.16)

• When no electronic cross-talk is assumed between sensors, the electronic gain matrix G is typically diagonal, i.e.,

G p = diag{g p }, (2.17)
meaning that each component is aected independently.

• A rotation matrix, as the Faraday rotation F i,p , is characterized by a given angle θ i,p and written as F i,p = cos(θ i,p ) -sin(θ i,p ) sin(θ i,p ) cos(θ i,p ) .

(2.18)

Let us note that the choice of coordinate system is important here since a rotation matrix becomes diagonal in a circular coordinate system [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF].

Jones matrices have their own specic variations w.r.t. time, frequency and direction. Usually, eects on the left hand side of the Jones chain in (2.13) are the same for all sources within the eld-of-view. Therefore, a matrix J i,p , for the specic propagation path i -p, can be decomposed into direction independent matrices Q p and direction dependent eects M i,p , leading to

J i,p = Q p M i,p K i,p (2.19)
where the scalar propagation delay K i,p clearly appears since it is the one fundamental eect to remain in an interferometer even when observations are completely uncorrupted. Thus, we clearly distinguish direction independent eects linked to the antennas and the design of the instrumentation system, from direction dependent perturbations as the beam pattern, Faraday rotation, parallactic angle, propagation through the tropospheric and ionospheric layers [START_REF] Noordam | The MeqTrees software system and its use for thirdgeneration calibration of radio interferometers[END_REF][START_REF] Rau | Advances in calibration and imaging techniques in radio interferometry[END_REF],

which are the most signicant perturbations to calibrate in advanced telescopes nowadays [START_REF] Bhatnagar | Correcting direction-dependent gains in the deconvolution of radio interferometric images[END_REF][START_REF] Lochner | Bayesian inference for radio observations[END_REF].

As regards time and frequency, scale variations dier as well according to the Jones matrix.

To give an example, the B p bandwidth matrix mainly depends on frequency while the electronic gain G p matrix is varying temporally essentially [START_REF] Noordam | The MeqTrees software system and its use for thirdgeneration calibration of radio interferometers[END_REF]. Let us mention that the ionospheric phase CHAPTER 2. PROBLEM SETUP AND STATE-OF-THE-ART delay is inversely proportional to frequency while Faraday rotation angle is inversely proportional to the square of frequency. These variations will be exploited in Chapter 4.

As it can be noticed from above, Jones matrices are a suitable representation to describe interferometric systems with dierent antennas and various propagation paths. Therefore, calibration amounts to estimate these specic matrices. Two options are possible in the estimation process: consider the combined eects in the global matrix J i,p and directly estimate its entries (this is the non-structured case [14]) or estimate unknowns in a parametric model and introduce structured matrices as in the Jones chain (2.14) (the so-called structured case [START_REF] Lochner | Bayesian inference for radio observations[END_REF]). Considering a parametric model in the structured scenario can lead to misspecications if the model does not correspond exactly to reality and estimation performance can easily deteriorate. However, such models introduce a few parameters to estimate and a smaller variance can be attained in the absence of misspecications. On the other hand, the non-structured case can be more robust as no particular specication is made about the model but far more unknowns need to be estimated, leading to possibly worse statistical behavior. In this work, we deal with both cases and depending on the operating environment, the user can consider either scenario.

In the next section, we expose state-of-the-art calibration methods to correct all amplitude and phase errors due to the environment and the instrumentation, described by the previously introduced Jones matrices.

Calibration methods

General strategies have been developed in radio astronomy to reach accurate calibration and thus, high dynamic-range images. Two classical approaches are presented hereafter.

Calibration strategies: external and self-calibration

• External calibration: Conventional calibration consists in estimating unknown gains while observing an isolated and bright calibration source whose structure, i.e., position and intensity, is given by tables. Source 3C286 is an example of external calibrator, used in the case of the VLA [START_REF] Leahy | VLA images of 23 extragalactic radio sources[END_REF]. In practice, calibration w.r.t. a reference source needs to be repeated frequently in order to know the evolution of distortions through time. Besides, the calibration source must be close enough to the unknown target (the source of interest) s.t. incurred disturbances are similar, otherwise, interpolation is necessary [START_REF] Yatawatta | On the interpolation of calibration solutions obtained in radio interferometry[END_REF]. Still, they must also be separated enough s.t. their contributions are clearly distinct in the data sets. Such method has some disadvantages: accurate a priori knowledge is required about properties of the calibration source and for antennas with wide eld-of-view, information is only provided in the direction of the calibration source.

• Self-calibration: Also called self-cal, this approach is an enhancement of external calibration.

Iteratively and alternatively, both the sky (parameters of interest) and dierent perturbations (nuisance parameters) are estimated, the two being mutually dependent. The self-cal procedure used to apply traditionally in regimes 1 and 2: there are no direction dependent eects but they may dier per antenna due to the instrumentation. This means that the eld-of-view is narrow and the same part of the ionosphere is observed. However, extension to the direction dependent calibration problem has also been addressed [START_REF] Van Der Tol | Self-calibration for the LOFAR radio astronomical array[END_REF].

The actual calibration and imaging strategy based on the self-cal and the LOFAR operation [START_REF] Jongerius | An end-to-end computing model for the Square Kilometre Array[END_REF][START_REF] Nijboer | LOFAR calibration[END_REF] is represented in Figure 2.6: the sky model is upgraded in major and minor cycles while calibration parameters are optimized by a designed algorithm in the calibration cycle, which is the subject of this document. More specically, thanks to a given sky model, we start from an To sum it up, the aim of self-cal was originally to estimate direction independent complex gains of antennas and noise parameters [START_REF] Friedlander | Direction nding in the presence of mutual coupling[END_REF][START_REF] Ng | Sensor-array calibration using a maximum-likelihood approach[END_REF]. But the problem has also been extended to multiple sources in the eld-of-view [START_REF] Fuhrmann | Estimation of sensor gain and phase[END_REF][START_REF] Wijnholds | A multisource calibration method for phased array radio telescopes[END_REF], where the assumed known source covariance matrix is equivalent to the source coherency matrix introduced in (2.10). Generalization to direction dependent perturbations, due to beam pattern or ionospheric eects, is also practicable [START_REF] Van Der Tol | Self-calibration for the LOFAR radio astronomical array[END_REF][START_REF] Noordam | Generalized self-calibration for LOFAR[END_REF].

In this specic case, some assumptions are needed about the structure of the ionospheric layer or geometry of the array.

Finally, improved variants of the self-cal procedure are numerous: (demixed) Peeling [START_REF] Noordam | LOFAR calibration challenges[END_REF][START_REF] Van Der Tol | Self-calibration for the LOFAR radio astronomical array[END_REF][START_REF] Jes | Direction dependent self calibration of large distributed sensor arrays[END_REF] is adapted for direction dependent gains and calibration is performed w.r.t. the brightest source in the sky which is then subtracted from the data and the procedure is repeated for the next brightest source. Such technique has notably been deployed for the VLA and the MWA [START_REF] Van Der Tol | Bayesian estimation for ionospheric calibration in radio astronomy[END_REF][START_REF] Mitchell | Real-time calibration of the Murchison Wideeld Array[END_REF]. We can also mention Pointing self-cal [START_REF] Bhatnagar | Solving for the antenna based pointing errors[END_REF], which corrects for errors in pointing directions, leading otherwise to distortions in the recovered image. In addition, the Field-Based Calibration (FBC) is suitable for ionospheric calibration in regime 3 and estimates phase gradients [10, 83 85] using interpolation and Zernike polynomials. Indeed, in regime 3, antennas are relatively close to each other so the ionospheric variation within the eld-of-view can be modeled as a linear gradient, using smooth polynomial model with low order. To generate such polynomials, Zernike or Karhunen-Loève transformations can be exploited. Another variant of Peeling is the Source Peeling and Atmospheric Modeling (SPAM) [START_REF] Intema | Ionospheric calibration of low frequency radio interferometric observations using the peeling scheme-i. method description and rst results[END_REF] which estimates ionospheric phase delays in regime 4. In this regime, antennas are well separated from one another and high order polynomials are more suitable to model the ionosphere.

CHAPTER 2. PROBLEM SETUP AND STATE-OF-THE-ART

In the next section, we give more details about array signal processing methods, which can be directly incorporated into the calibration cycle of Figure 2.6.

LS calibration

One of the most popular and statistically ecient method in signal processing is the ML estimator. However, it may suer from heavy computational cost and it is not solvable in closed-

form. An alternative option is to consider Least Squares (LS) based algorithms [START_REF]Multisource self-calibration for sensor arrays[END_REF][START_REF] Van Der Veen | Signal processing tools for radio astronomy[END_REF] which are asymptotically ecient for large number of samples under Gaussian noise [START_REF] Ottersten | Covariance matching estimation techniques for array signal processing applications[END_REF]. Besides, they benet from reduced complexity w.r.t. the original ML and consist in minimizing the sum of squares of errors between data points and the predicted model. A possible nonlinear optimization technique to solve such problems is the gradient-based Levenberg-Marquardt (LM) algorithm [START_REF] Madsen | Methods for non-linear least squares problems[END_REF], exposed briey in the following.

LM algorithm

Using noisy observations from (2.12), its vectorization leads to

vec(V pq ) = v pq = D i=1 J * i,q (θ) ⊗ J i,p (θ)vec(C i ) + n pq (2.20)
with n pq = vec(N pq ). Jones matrices are variable with time and frequency but as a rst step, we are considering a suciently ne interval so that variations are neglected. Time t and frequency f dependence are thus not specied for now.

Considering only cross-correlations for all antenna pairs [14, 89], one obtains the following 4B × 1 vector, with B = M (M -1)/2 the total number of baselines,

y = [v 12 , v 13 , . . . , v (M -1)M ] T = D i=1 s i (θ) + n (2.21)
where s i incorporates the contribution of known coherency of the i-th source and the Jones matrices for all baselines, i.e., s i (θ) = s T i,12 (θ), s T i,13 (θ), . . . , s T i,(M -1)M (θ)

T and s i,pq (θ) = J * i,q (θ) ⊗ J i,p (θ)vec(C i ).

The corresponding expression of the LS cost function is given by

θ = argmin θ ||y -s(θ)|| 2 2 (2.22) in which s(θ) = D i=1 s i (θ). If we note ζ(θ) = ||y -D i=1 s i (θ)|| 2
2 the considered cost function, the (t +1)-th iteration of the LM algorithm reads

θ t+1 = θ t -∇ θ ∇ T θ ζ(θ) + λdiag{∇ θ ∇ T θ ζ(θ)} -1 ∇ θ ζ(θ)| θ t . (2.23)
in which ∇ θ is the gradient with respect to θ and ∇ θ ∇ T θ is the Hessian matrix. Numeri- cal implementation of the LM is exposed in [START_REF] Madsen | Methods for non-linear least squares problems[END_REF][START_REF] Gavin | The Levenberg-Marquardt method for nonlinear least squares curve-tting problems[END_REF]. Let us note that a possible adaptation of the LM considers an identity matrix instead of the diagonal of the Hessian matrix, noted diag{∇ θ ∇ T θ ζ(θ)} here. Any adapted positive-denite diagonal matrix can be used.

The LM algorithm [START_REF] Levenberg | A method for the solution of certain nonlinear problems in least squares[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] is based on a damped Gauss-Newton method, i.e., a combination between two minimization techniques: the gradient descent and the Gauss-Newton, depending 2.4. CALIBRATION METHODS

21

on how far the estimated parameters are from the optimal values. Choice of the damping parameter λ and initializations impact convergence and reaching the global optimum is not always possible. Due to computation of a Jacobian and inversion of large-dimensional matrices in (2.23), convergence can be extremely slow. The computational cost gets even more signicant if the number of parameters to estimate is huge, as it is the case for very large radio telescopes.

Speeding up the procedure is possible using Graphics Processing Unit (GPU) [START_REF] Yatawatta | GPU accelerated nonlinear optimization in radio interferometric calibration[END_REF] or simplifying the problem with variants of the LS as EM or the Space-Alternating Generalized Expectation-Maximization (SAGE) algorithms, presented in section 2.4.3.

In the following section, we aim to present the (Weighted) Alternating Least-Squares ((W)ALS) algorithm which is asymptotically unbiased and ecient. It leads to closed-form expressions for some of the estimates and fullls the convergence properties.

(W)ALS algorithm

In this subsection, we present array signal processing methods, based on the so-called covariance model. Such model is slightly dierent from the Jones model in (2.10) and was initially used to design the ALS and its variants. Nevertheless, the same approach can be applied to the Jones-based model.

Let us consider M antennas in an array with impinging D signal sources. Each antenna p with p ∈ {1, . . . , M } is described by its position l p = [x p , y p , z p ] T which is known. Likewise, each source i with i ∈ {1, . . . , D} is dened by its specic spatial location

d i = [l i , m i , n i ] T s.t. n i = 1 -l 2 i -m 2 i [60]
. The M × D steering matrix is thus given by where s and n refer to the source and noise contribution, respectively. The noise is usually assumed to follow a zero-mean white Gaussian distribution. From a signal processing perspective, signal sources can be modeled in two dierent ways: either they are assumed deterministic or stochastic, i.e., random and following a given a priori distribution.

A = exp -j 2π λ LD (2.
In our specic application, interferometers compute correlations so we consider covariance matching estimation techniques. In the ideal case, when no perturbation occurs, except for the noise, the resulting covariance matrix reads [START_REF] Wijnholds | Signal processing challenges for radio astronomical arrays[END_REF] R = E r(t)r

H (t) = ASA H + N (2.26)
where S is the source covariance matrix (known as the coherency/brightness matrix in the Jonesbased model) with source powers along the diagonal and N is the noise covariance matrix which is diagonal in classical cases, with variable or not noise powers along the diagonal. Thus, the set of observations consists in covariance matrices which are functions of the source structure, the environmental and instrumental errors, and the receiver noise. The ALS approach is a suitable technique for solving LS covariance model tting as follows

θ = argmin θ || R -R(θ)|| 2 F (2.27)
where R is the sample covariance matrix, i.e., an estimate of the data covariance matrix, given by • Regime 1: In the general case, calibration requires to estimate direction independent elec- tronic gains G = diag{g} s.t.

T -1 YY H with Y = [r(t
R = GASA H G H + N.
(2.28)

These gains are the same for all sources in the eld-of-view which is narrow. Consequently, in regime 1, we only need to observe a calibration source whose intensity and position are given.

If we are pointing towards this specic source, the matrix A reduces to a vector a and S to the scalar source power σ 2 s which are both known from tables, leading to the simplied problem

R = gg H + N (2.29)
where g = Gm with m = aσ s , known as a rank-1 factor analysis model from which g and N can be solved [START_REF] Boonstra | Gain calibration methods for radio telescope arrays[END_REF].

• Regime 2: The ionospheric perturbations are still the same for the entire eld-of-view due to its narrowness but each individual antenna is assigned a unique perturbation term, leading to a similar calibration procedure as in Regime 1.

• Regime 3: Direction dependent eects and beam patterns of antennas need to be taken into account. To this end, we introduce direction dependent gains γ s.t.

R = GAΓSΓ H A H G H + N (2.30)
where Γ = diag{γ}. For the calibration sources, S is known and we usually consider Σ = ΓSΓ H as a global unknown diagonal matrix [START_REF] Wijnholds | A multisource calibration method for phased array radio telescopes[END_REF].

• Regime 4: A gain needs to be estimated for each source and each antenna, since many unknown sources are visible in the eld-of-view. The data measurements are given by [START_REF] Van Der Tol | Self-calibration for the LOFAR radio astronomical array[END_REF]

R = (G A)S(G A) H + N (2.31)
where G is no longer diagonal, but full and considers all possible perturbation eects whether they are direction dependent or not. Such calibration problem requires further assumptions in order to be identiable.

Advanced interferometric arrays are described by regimes 3 and 4 in which directional eects are prevailing. Therefore, we will focus much of our attention on these two regimes:

• In regime 3, we need to estimate direction independent instrumental gains G, the source powers Σ = diag{σ} and the noise powers N = diag{σ n } as follows [START_REF]Multisource self-calibration for sensor arrays[END_REF][START_REF] Wijnholds | Fish-eye observing with phased array radio telescopes[END_REF]

ĝ = argmin g || R -GAΣA H G H -N|| 2 F (2.32) σ = argmin σ || R -GAΣA H G H -N|| 2 F (2.33) σn = argmin σn || R -GAΣA H G H -N|| 2 F .
(2.34) Equation (2.32) requires a rank-1 approximation in order to get an estimate while (2.33) and

(2.34) lead to closed-form expressions. In fact, the gain estimation step can be addressed in various ways: iteratively with a Gauss-Newton method (GNLS) [START_REF] Gill | Practical optimization[END_REF], logarithmically with the Logarithmic Least Squares (LOGLS) [START_REF] Boonstra | Gain calibration methods for radio telescope arrays[END_REF] or even using Column Ratios (COLR) [START_REF]Multisource self-calibration for sensor arrays[END_REF][START_REF] Wijnholds | Fish-eye observing with phased array radio telescopes[END_REF].

Powers and positions of the calibration sources, as well as antenna locations are typically assumed known. Due to the possible large number of unknowns to estimate, constraints can be imposed on a reference source with xed power, phase and position in order to avoid any identiability issues. However, due to ionospheric phase shifts, exact positions of calibration sources may still be known inaccurately. To address this problem, a Direction Of Arrival (DOA) estimation step can be added in the iterative procedure. These inaccurate location parameters appear in A whose parametric model is available, and are estimated with MUSIC or Weighted Subspace Fitting (WSF) algorithms [START_REF] Viberg | Detection and estimation in sensor array processing using Weighted Subspace Fitting[END_REF], which are asymptotically statistically ecient and make use of signal and noise subspaces [START_REF] Haardt | Subspace methods and exploitation of special array structures[END_REF]. Incorporating the DOA estimation in the iterative procedure leads to an extension of ALS, named xALS [START_REF]Multisource self-calibration for sensor arrays[END_REF][START_REF] Wijnholds | Fish-eye observing with phased array radio telescopes[END_REF].

• In regime 4, the number of unknowns is tremendous. To be able to solve the problem and remove identiability issues, physical constraints must be added. These constraints can be related to the structured array into multiple compact subarrays (regime 3) or the specic variation of perturbations w.r.t. time and frequency, with a dierence between instrumental gains and ionospheric eects. Imposing constraints enables to reduce the number of unknown parameters to estimate so that the Fisher Information Matrix (FIM), whose inverse leads to the CRB, is invertable with no singularities.

If we exploit time t and frequency f , the model (2.31) can be written as [START_REF] Van Der Tol | Self-calibration for the LOFAR radio astronomical array[END_REF] R

[f ] t = (G [f ] t A [f ] t )S [f ] t (G [f ] t A [f ] t ) H + N [f ] t (2.35)
and the calibration problem is considered for each time interval and frequency band, as follows

θ[f] t = argmin θ [f ] t || R[f] t -R [f ] t (θ [f ] t )|| 2 F s.t. C {θ [f ] t } t,f = 0 (2.36)
where C stands for the considered constraints w.r.t. time and frequency.

In the following, we describe the Statistically EFcient and Fast Calibration (StEFCal) method, adapted for antenna-based gains calibration [START_REF] Salvini | Fast gain calibration in radio astronomy using alternating direction implicit methods: Analysis and applications[END_REF][START_REF]StEFCal an alternating direction implicit method for fast full polarization array calibration[END_REF], which provides computational advantage over the previously introduced algorithms derived from (weighted) LS cost function.

CHAPTER 2. PROBLEM SETUP AND STATE-OF-THE-ART

StEFCal algorithm

Introduced in [START_REF] Salvini | Fast gain calibration in radio astronomy using alternating direction implicit methods: Analysis and applications[END_REF], this Alternating Direction Implicit (ADI) approach proposes to reduce the numerical complexity of typical LM solver which is restraining for implementation in modern radio interferometers. Its convergence to optimal solution is guaranteed and it reveals to be more ecient than traditional antenna-based gain methods.

In this method, the gain of each receive path is estimated, assuming that gains of all other antennas are already calibrated. The iterative procedure is then repeated successively for each antenna. Such methods have already been considered in [START_REF]StEFCal an alternating direction implicit method for fast full polarization array calibration[END_REF][START_REF] Mitchell | Real-time calibration of the Murchison Wideeld Array[END_REF][START_REF] Repetti | Non-convex optimization for selfcalibration of direction-dependent eects in radio interferometric imaging[END_REF]. To improve convergence, robustness and reduce the inuence of noise on the estimates, numerical simulations have shown that averaging even and odd iterations accelerates the convergence. Powers and positions of sources being known, unknown parameters to estimate reduce to the complex gains describing the elements of the array and receiver path noise powers, leading to the following calibration problem {ĝ, σn } = argmin g,σn

||Ω H ( R -GASA H G H -N)Ω|| 2 F (2.37)
where Ω is a possible weighting term.

Autocorrelations being essentially dominated by noise, estimation of gains is performed using cross-correlations only. Thus, to simplify the problem during the gain calibration step, the noise N is ignored and the diagonals of R and M are set to 0. We note M the available model of the observed eld, which is approximately equal to ASA H and only includes the undistorted brightest calibration sources. The estimation problem is then described by

ĝ = argmin g || R -GMG H || 2 F . (2.38)
The underlying approach consists in estimating g H for xed g, and vice versa. Considering decoupled complex gains, the objective function in (2.38) can be written as

|| R -GMG H || 2 F = M p=1 ||[ R] :,p -[GMG H ] :,p || 2 F .
(2.39)

Therefore, the problem is divided into independent linear LS problems and each sensor gain is updated successively for p ∈ {1, . . . , M } during one given iteration.

Statistical estimation performance is similar whether we use directly the LM as an optimization technique or the StEFCal procedure. Still, this last method results in lower complexity since

O(M 2 ) is required instead of O(M 3
). Variants of StEFCal are exposed in [START_REF] Salvini | Fast gain calibration in radio astronomy using alternating direction implicit methods: Analysis and applications[END_REF]. For instance, to reduce the eects of outliers, a ponderation term can be incorporated leading to the Iterative Reweighted Least Squares (IRLS).

The StEFCal method can also be directly used as the gain estimation step of the ALS procedure in (2.32) so as to improve the implementation, the computational cost and the speed of convergence. Indeed, gain calibration reveals to be almost 19 times faster with StEFCal than with the WALS method described in [START_REF]Multisource self-calibration for sensor arrays[END_REF]. Still, StEFCal was developed for estimation of direction independent gains, which is limited.

Considering all the specic introduced perturbations and their corresponding statistics enables to exploit all available information. In the following section, the Maximum A Posteriori (MAP) estimator is notably based on a priori knowledge.

Bayesian algorithm

Bayesian approaches are alternatives to ML-based techniques and one example is the MAP estimator [START_REF] Van Der Tol | Bayesian estimation for ionospheric calibration in radio astronomy[END_REF][START_REF] Van Der Tol | Ionospheric calibration for the LOFAR radio telescope[END_REF], in which prior information is taken into account. As an example, the a priori model for the ionosphere should be accurate enough to reect the uctuations on a very small scale but not too complex so that the estimation algorithm can be deducted at low cost. In this sense, a statistical model can be introduced to describe the phase uctuations within the ionosphere [START_REF] Van Der Tol | Bayesian estimation for ionospheric calibration in radio astronomy[END_REF], notably high-order statistics in order to study local dierences. Kolmogorov theory on turbulent environment is particularly well-adapted. Let us note that introducing many parameters to describe the model makes it more accurate and realistic, leading to small bias. On the other hand, with less parameters to estimate, we can hope to achieve a better variance, i.e., a lower estimation error.

One of the most popular estimators in estimation theory is the Minimum Mean Squared Error (MMSE), given by

θ = argmin θ E || θ -θ|| 2 2 , (2.40)
i.e., the minimization of the Mean Square Error (MSE). The corresponding solution is θ = E θ|y , so posterior mean of the parameters needs to be computed. According to Bayes's law, this distribution can be written as p(θ|y) = p(y|θ)p(θ) p(y|θ)p(θ)dθ .

(2.41)

The MMSE is the optimal estimator in the LS sense but remains dicult to implement due to multidimensional integral and expected value from (2.41) and (2.40).

A good alternative is the MAP estimator which consists in maximizing the a posteriori distribution, i.e., (2.42)

When the noise follows a Gaussian distribution, this amounts to solve a non-linear LS problem.

In our interferometric calibration problem, the MAP estimator leads to [START_REF] Wijnholds | Calibration challenges for future radio telescopes[END_REF][START_REF] Van Der Tol | Bayesian estimation for ionospheric calibration in radio astronomy[END_REF] θ = argmax

θ ||Ψ -1 2 (y -s(θ)) || 2 2 + ||C -1 2 θ θ|| 2 2 (2.43)
where Ψ is the covariance of the noise vector n and we note C θ the covariance for θ. In the case of a pure Kolmogorov turbulent process, C θ can be modeled by a power law with unknown parameters, thus, a priori knowledge is exploited and the additional term in (2.43) acts as a penalty function. Let us note that from a Bayesian perspective, the introduced hyperparameters in the model for C θ must be estimated as they describe its statistical structure. However, the main issue is to obtain accurate models for the unknown perturbations, which are not always available in practice.

Variants to LS calibration

LS calibration represents the traditional estimator in case of additive white Gaussian noise.

Still, there are some well-known disadvantages: speed of convergence is slow, convergence to global optimum is not guaranteed and the computational cost is signicant. To improve both complexity and accuracy, alternative methods have been introduced.

EM algorithm

The EM algorithm [START_REF] Feder | Parameter estimation of superimposed signals using the EM algorithm[END_REF][START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is an iterative method which provides ML estimates asymptotically when observations are viewed as incomplete data, i.e., some data points are unknown because they are missing or hidden. Compared to classical optimization algorithms used to solve non-linear LS problems, as the LM, the EM enables to reduce the computational cost and speed up the convergence.

It consists of two consecutive steps: an expectation (E) and a maximization (M) steps, and requires to dene the so-called complete data w, extracted from the output data y. In order to specify the complete data, in our context, we assume well separated sources with unique propagation paths, from one source to one antenna, so that we can assign one unknown parameter vector θ i to each source for i ∈ {1, . . . , D}. Thus, a partitioning per source can be conducted.

The case of sources sharing some parameters, due for instance to small angular deviation, is addressed in [START_REF] Kazemi | Radio interferometric calibration using the SAGE algorithm[END_REF].

Instead of considering the complete data model with contributions from all directions, we reduce the problem to multiple sub-problems, leading to the following complete data vector for the i-th source [14, 89]

w i = s i (θ i ) + n i (2.44)
where the unknown parameter vector is decomposed as θ

= [θ T 1 , • • • , θ T D ]
T . The noise is also assumed to be decomposable into source contributions, i.e., n = D i=1 n i and y = D i=1 w i . To rene the description of the model, it is possible to associate more or less noise to each source, depending on its brightness, thanks to a weighting term.

The two steps of the EM algorithm are summarized as follows

• E step: in our case, it amounts to compute the expectation value of complete data condi- tionally to the observed data and the unknown parameter vector, i.e., ŵ = E w|y; θ .

(2.45)

In the Gaussian case, a closed-form expression can be attained [START_REF] Anderson | An introduction to multivariate statistical analysis[END_REF].

• M step: it consists in minimizing a cost function, for one given source i and repeat the estimation for i ∈ {1, . . . , D} so the procedure is parallelisable. The cost function under study is thus given by ζ i (θ

i ) = ||w i -s i (θ i )|| 2 F , instead of ζ(θ) = ||y -D i=1 s i (θ)|| 2
F , which enables to reduce the complexity of the problem. To perform numerical minimization, the LM algorithm can still be used and computation is performed with respect to θ i instead of θ.

SAGE algorithm

The SAGE algorithm [START_REF] Fessler | Space-alternating generalized expectation-maximization algorithm[END_REF] is an alternative method to the EM. It can be applied for directiondependent calibration of radio interferometers, leading to the SAGECal algorithm, whose implementation time and convergence speed are better than the EM. Compared to the latter, a dierent denition of the complete data is considered and the noise is assigned dierently.

Specically, the hidden data are selected as [14]

w i = s i (θ i ) + n (2.46) s.t. y = w i + D k=1 k =i s k .
The two steps of the SAGE are quite similar to those of the EM: they consist in computing a conditional mean followed by a minimization process.

The SAGE algorithm is faster and leads to better calibration accuracy than the EM, but they are both superior in numerical complexity compared to direct application of the LM on the LS problem. To resolve issues on a more timely basis, it is therefore relevant to reduce the dimension of the unknown vector so as to invert matrices with lower dimensions. To ensure convergence towards a global optimum, probabilistic techniques as the Simulated Annealing (SA) [START_REF] Kirkpatrick | Optimization by simulated annealing: Quantitative studies[END_REF] are possible but on the other hand, they unfortunately reduce the convergence speed.

The SAGE algorithm is particularly appropriate in some pathological cases when sources are not easily dierentiable from the background noise or when they are close to each other [START_REF] Chung | Array self-calibration using SAGE algorithm[END_REF] or even when positions of sensors are misspecied. Generally, classical calibration sources are viewed as bright and compact point sources with no broad or diuse structure. But they may not be bright enough. To remedy this, we can introduce clustered sources [START_REF] Kazemi | Clustered calibration: an improvement to radio interferometric direction-dependent self-calibration[END_REF] where each clustering acts as a single bright source. More information is available as more sources are considered all together and not individually. To this end, sources of the same cluster must be physically close to one another with a weighting to take into account their respective intensity and the perturbations must be similar due to small angular deviation. Such technique simplies the problem since we reduce the number of directions to calibrate and we consider the information associated to several sources all together.

All the previously mentioned methods are applied assuming a Gaussian noise model for n pq in (2.20) in the context of radio interferometry. Still, this is not always realistic and robustness needs to be addressed.

ECME algorithm

Due to the presence of outliers which corrupt the observations, the noise cannot simply be considered Gaussian [START_REF] Kazemi | Robust radio interferometric calibration using the tdistribution[END_REF]. These are due to radio frequency interferences which are not always perfectly ltered [START_REF] Raza | Spatial ltering of RF interference in radio astronomy[END_REF][START_REF] Leshem | Multichannel interference mitigation techniques in radio astronomy[END_REF][START_REF] Leshem | Radio-astronomical imaging in the presence of strong radio interference[END_REF], to errors in the sky model due to weak unknown sources in the background [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF], to the temporary appearance of particularly bright sources in the sidelobes generating punctual interferences [START_REF] Boonstra | Radio frequency interference mitigation in radio astronomy[END_REF].

To investigate robustness [START_REF] Zoubir | Robust estimation in signal processing: A tutorial-style treatment of fundamental concepts[END_REF], the noise can be modeled as a Student's t distribution [START_REF] Lange | Robust statistical modeling using the t distribution[END_REF].

Such distribution is a suitable heavy-tailed statistical model for data sets aected by the presence of outliers. To provide relevant ML estimates with a Student's t noise model, the Expectation-Conditional Maximization Either (ECME) algorithm [START_REF] Liu | ML estimation of the t distribution using EM and its extensions[END_REF], an extension of the EM, is considered in [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF][START_REF] Kazemi | Robust radio interferometric calibration using the tdistribution[END_REF] due to its adequate rate of convergence. Using (2.21), we distinguish the contribution of D known bright sources from the contribution of noise and D unknown weak sources. The latter act as outliers and can be incorporated in a total non-Gaussian noise n, s.t.,

y = D s=1 s i (θ) + n (2.47) in which n = D i =1 s i (θ )
ST (t; µ, σ 2 n , γ) = Γ( ν+1 2 ) √ πνΓ( ν 2 )σ n 1 + 1 ν t -µ σ n 2 -ν+1 2 .
(2.48)

As the degree of freedom approaches innity, the Student's t becomes a Gaussian distribution while it is heavy-tailed for low values of ν [START_REF] Bartkowiak | Should normal distribution be normal? The Student's t alternative[END_REF]. If n is assumed to follow the xed distribution in (2.48) with σ n = 1, µ = s(θ) and independent and identically distributed (i.i.d.) entries, the corresponding likelihood expression is a product between N components where N is the number of real data points in y. In this specic case, maximizing the resulting likelihood does not provide the ML estimate as the noise is not Gaussian and ν needs to be estimated as well.

In the ECME, auxiliary variables are introduced, such as weights w and a scalar λ. Two steps are then considered:

• The E step consists in computing the following conditional expression ŵ = E w|y; θ, ν

for each entry of the hidden vector w. The scalar λ is also updated as function of the weights w.

• The goal of the maximization step is to provide an estimate of ν by maximizing the loglikelihood expression w.r.t. ν and equate the result to zero. But the most important part is to obtain an estimate for θ. To this end, the following minimization process is considered

θ = argmin θ N k=1 [w] k [y] k -[s(θ)] k 2 (2.50)
with a weighted LS cost function. To solve (2.50), a robust version of the LM can be considered [START_REF] Kazemi | Robust radio interferometric calibration using the tdistribution[END_REF].

The robustness issue has been tackled but all the previously introduced algorithms consider one frequency bin at a time during calibration, which is not an optimal way to process the data. Therefore, in the following, we focus on a multi-frequency scenario considering the full observation bandwidth as a whole.

ADMM-based algorithm

To exploit frequency dependence of Jones matrices, it is possible to introduce polynomials, with a given order, assuming smooth variation across frequency. In a multi-frequency scenario, the computing load is particularly heavy due to large volume of data. To relieve this burden, data parallelism can be brought to bear by using a network of multiple agents among which data is distributed across frequency [START_REF] Yatawatta | Distributed radio interferometric calibration[END_REF]. Therefore, data are not centralized but relocated in space and each agent has access to some specic data, at a given frequency. Calibration can then be reformulated as a consensus distributed optimization problem [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. As a rst step, we only consider frequency dependence here, in a given time interval.

For each frequency f , Jones matrices are modeled as follows to enforce smoothness

J [f ] = P [f ] Z (2.51)

CALIBRATION METHODS

29

where P [f ] collects all the frequency-dependent polynomial terms and Z refers to the hidden variables which need to be estimated.

Instead of estimating directly Z from all data available, distributed calibration consists in performing calibration locally by each agent, for a given subband. Thus, calibration amounts to solve the following constrained problem

{ Ĵ[f 1 ] , Ĵ[f 2 ] , • • • , Ẑ} = argmin J [f 1 ] ,J [f 2 ] ,••• ,Z f ψ [f ] (J [f ] ) s.t. J [f ] = P [f ] Z (2.52) in which ψ [f ] (•) is the following LS cost function ψ [f ] (θ [f ] ) = pq ||V [f ] pq - D i=1 J [f ] i,p (θ [f ] )C [f ] i J [f ] H i,q (θ [f ] )|| 2 F .
(2.53)

To solve (2.52), we introduce unknown Lagrange multipliers Y [f ] and a xed regularization term ρ leading to the following Lagrangian

L(J [f 1 ] , J [f 2 ] , • • • , Z, Y [f 1 ] , Y [f 2 ] , • • • ) = f L [f ] (J [f ] , Z, Y [f ] ) (2.54) in which L [f ] (J [f ] , Z, Y [f ] ) = ψ [f ] (J [f ] ) + ||Y [f ] H J [f ] -P [f ] Z || F + ρ 2 ||J [f ] -P [f ] Z|| 2 F . (2.55)
A possible method to solve this problem is the Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF][START_REF] Li | Alternating Direction Method of Multipliers for separable convex optimization of real functions in complex variables[END_REF][START_REF] Jiang | Alternating direction method of multipliers for real and complex polynomial optimization models[END_REF]. Instead of solving a global problem as argmin f ψ [f ] (•) at the fusion center which collects data from all agents, calibration is performed at each agent locally, for one given frequency. Therefore, frequencies are processed independently but simultaneously and communication between computational agents is ensured thanks to the global fusion center. The study was originally introduced for one single source but it can easily be extended to the case of D calibration sources by using parameter partitioning with SAGE algorithm [START_REF] Kazemi | Radio interferometric calibration using the SAGE algorithm[END_REF].

The three iterative steps of the ADMM are the following

Ĵ[f] t+1 = argmin J [f ] L [f ] θ [f ] , Ẑ t , Ŷ[f] t (2.56)
performed locally by each agent

Ẑ t+1 = argmin Z f L [f ] Ĵ[f] t+1 , Z, Ŷ[f] t (2.57)
performed at the fusion center

Ŷ[f] t+1 = Ŷ[f] t + ρ Ĵ[f] t+1 -P [f ] Ẑ t+1 (2.58)
performed locally by each agent where t is the iteration counter. Estimation in (2.56) can be done iteratively with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) [START_REF] Nocedal | Numerical optimization[END_REF] or the Riemannian Trust-Region (RTR) [START_REF] Absil | Trust-region methods on Riemannian manifolds[END_REF] algorithms.

CHAPTER 2. PROBLEM SETUP AND STATE-OF-THE-ART

Each Jones estimate Ĵ[f] obtained by each agent is transfered to the processing center which then estimates Z globally. Indeed, estimation of Z requires information across all frequencies and enforces consensus among all agents. Solving (2.57) leads to a closed-form expression. Finally, update of Y [f ] is performed in the last step by each corresponding agent. Let us note that similar techniques have been studied for image synthesis [START_REF] Carrillo | PURIFY: a new approach to radiointerferometric imaging[END_REF][START_REF] Ferrari | Distributed image reconstruction for very large arrays in radio astronomy[END_REF].

The advantages of distributed optimization in the network include the reduced computational cost as calculus are distributed among several agents. To go even further, it would be meaningful to take into account temporal variation of Jones matrices or treat extreme cases, such as when the number of frequencies is very high compared to the number of available agents. But all these studies are still under investigation.

In the last section of this chapter, we address some notions as regards the imaging process, without going into details as this is not the central core of the manuscript.

Imaging in radio astronomy

Measured correlations for baselines (u, v) and true brightness image I are related by the Van Cittert-Zernike theorem in (2.5). Producing an image, namely map making, consists in inverting such expression and leads to the dirty image as in (2.7) due to the nite number of observations.

The dirty image can also be written as a convolution between the desired image I and the Point Spread Function (PSF) B, also called the dirty beam, which is the inverse Fourier transform of the sampling function S, i.e.,

I D = I * B (2.59)
and B(l, m) = S(u, v)e 2πj(ul+vm) dudv.

(2.60)

In order to recover the nal restored image I from I D , with high-resolution, a priori knowledge about B and deconvolution techniques are required. Thus, the deconvolution process is essentially a Fourier inversion problem and consists in removing the instrumental response. Let us note that autocorrelations are usually not used in the image formation, in order to reduce the inuence of noise. Besides, forming an image with high dynamic range becomes all the more dicult to perform with the new generation of radio interferometers, the new challenges at stake and the signicant perturbation eects which need to be calibrated and removed.

Many deconvolution algorithms are based on sequential source removal and the most popular one is CLEAN [START_REF] Högbom | Aperture synthesis with a non-regular distribution of interferometer baselines[END_REF]: from a mostly empty sky model, composed of discrete point sources, it sequentially estimates the brightest source and subtracts its partial contribution from the dirty image using a so-called loop gain, until reaching a noise-like image where only the weakest possible sources remain. Iteratively, each source position and intensity are estimated and added to a point source list which will nally result in the global reconstructed image, convolved with an ideal reconstruction beam. Other alternative methods are possible, such as the Maximum Entropy Method (MEM) [START_REF] Cornwell | A simple maximum entropy deconvolution algorithm[END_REF], Clark CLEAN and Cotton-Schwab (CS) which accelerate and reduce the computational load of CLEAN [START_REF] Levanda | Synthetic aperture radio telescopes[END_REF], the parametric estimation based algorithm Least Squares Minimum Variance Imaging (LS-MVI) [START_REF] Ben-David | Parametric high resolution techniques for radio astronomical imaging[END_REF] or even sparse reconstruction techniques [START_REF] Carrillo | PURIFY: a new approach to radiointerferometric imaging[END_REF]121123]. Besides, since the spectral dimension can no longer be omitted, three-dimensional reconstruction approaches are increasingly becoming prevalent [START_REF] Rau | A multi-scale multi-frequency deconvolution algorithm for synthesis imaging in radio interferometry[END_REF]. As mentioned before, the traditional CLEAN is based on a point source model but extension with shapelets or wavelets is a possibility.

Lastly, let us note that two softwares have been implemented for calibration of the LOFAR [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF][START_REF] Nijboer | LOFAR calibration[END_REF] with direction dependent eects [START_REF] Bhatnagar | Calibration and imaging challenges at low radio frequencies: An overview of the state of the art, The Low Frequency Radio Universe[END_REF]: the MeqTrees software [START_REF] Noordam | The MeqTrees software system and its use for thirdgeneration calibration of radio interferometers[END_REF] which is able to simulate and calibrate radio astronomical measurements, and the Blackboard system [START_REF] Loose | LOFAR Self-Calibration using a Blackboard software architecture[END_REF][START_REF] Pandey | Calibrating LOFAR using the Black Board Selfcal System[END_REF], named Black-Board Self-cal (BBS). MeqTrees is particularly exible to specify the dierent perturbation eects to introduce and will be presented in Chapter 5 for realistic data sets.

Furthermore, traditional calibration and imaging were originally implemented in the Astronomical Image Processing System (AIPS) software package, which became AIPS++ [START_REF] Sault | The Hamaker-Bregman-Sault measurement equation[END_REF] and is now referred to as CASA [START_REF] Mcmullin | CASA architecture and applications[END_REF][START_REF] Myers | CASA synthesis & single dish reduction cookbook[END_REF]. Awimager is an example of imager, on CASA. Contrary to non directional eects which are usually corrected separately from the imaging process, estimation of errors due to direction dependency is usually combined with the deconvolution step [START_REF] Rau | Advances in calibration and imaging techniques in radio interferometry[END_REF][START_REF] Bhatnagar | Wide-eld wide-band interferometric imaging: The WB A-projection and hybrid algorithms[END_REF][START_REF] Bhatnagar | Correcting direction-dependent gains in the deconvolution of radio interferometric images[END_REF].

Therefore, reconstructed images will also be provided in this work to emphasize the benet of the designed direction dependent calibration techniques.

Conclusion

Most calibration algorithms are LS based, which operate iteratively for one given frequency channel and reveal to be asymptotically ecient under a Gaussian noise model. In section 2.4.3.3, the Student's t distribution was specically considered to include the contribution of interferences and weak background sources in radio interferometric data for non-Gaussian environment. In addition, the multi-frequency case was discussed in section 2.4.3.4 where the calibration task was reformulated as a consensus optimization problem. In the following chapters, we propose robust calibration algorithms based on a broader class of distributions so as to avoid model misspecications and achieve more robustness. An extension to the multi-frequency scenario is also presented and realistic simulations are performed with the exible software MeqTrees in order to reconstruct residual images.

Chapter 3

Robust calibration method in non-Gaussian environment

Motivations for robust calibration

Most calibration approaches are LS based [START_REF]Multisource self-calibration for sensor arrays[END_REF][START_REF] Boonstra | Gain calibration methods for radio telescope arrays[END_REF][START_REF] Salvini | Fast gain calibration in radio astronomy using alternating direction implicit methods: Analysis and applications[END_REF][START_REF] Wijnholds | Fish-eye observing with phased array radio telescopes[END_REF] and as mentioned in section 2.4.3, alternative algorithms as the EM [START_REF] Moon | The Expectation-Maximization algorithm[END_REF][START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF] or the SAGE [START_REF] Fessler | Space-alternating generalized expectation-maximization algorithm[END_REF] have been proposed to enhance the convergence rate. Still, the Gaussian noise assumption is usually considered in radio astronomy [14, 89], which is not always realistic and can lead to poor estimation performance. As explained in section 2.4.3.3, the presence of disruptive outliers has multiple causes: man-made radio signals generate frequency interferences which can be removed by agging [START_REF] Oringa | Postcorrelation radio frequency interference classication methods[END_REF] but are not always perfectly ltered in practice. The solar activity or the punctual appearance of strong sources in the sidelobes constitute other sources of interference which can randomly generate outliers. Furthermore, radio interferometric data include the contribution of bright calibration sources but also of multiple unknown weak sources in the background leading to incomplete sky models with ux loss and spurious sources [START_REF] Grobler | Calibration artefacts in radio interferometry I. Ghost sources in Westerbork Synthesis Radio Telescope data[END_REF][START_REF] Martí-Vidal | Spurious source generation in mapping from noisy phaseself-calibrated data[END_REF]. Thus, the data model for baseline (p, q) can be written as

V pq = D i=1 J i,p (θ)C i J H i,q (θ) + N pq (3.1)
where the noise component N pq includes the ambient noise N pq , which is assumed Gaussian, but also the presence of all previously mentioned sources of interferences leading to outliers.

Considering all available baselines leads to (2.47) with n and n, respectively, the total Gaussian and non-Gaussian noise components. The aim of robust calibration is thus to estimate θ under non-Gaussian noise environment n [START_REF] Meriaux | Robust-COMET for covariance estimation in convex structures: Algorithm and statistical properties[END_REF][START_REF]Ecient estimation of scatter matrix with convex structure under t-distribution[END_REF]. Let us note that during the imaging process, the D weak unknown sources, hidden in n, become the key elements of interest we wish to recover.

To robustify the calibration scheme, a distribution class, dierent from the Gaussian one, needs to be considered to model the noise contribution n. To the best of our knowledge, the only robust calibration procedure was proposed in [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF] and exposed in section 2.4.3.3. In this case, the noise part is specically modeled as a Student's t distribution with i.i.d. entries, which is not optimal and leads to model misspecications as the true distribution of n is not known in practice. In this chapter, we aim to show the limitations of the classical Gaussian noise assumption and propose an estimator based on a broader class of distributions. The noise vector n pq in (2.20) is traditionally assumed to follow a Gaussian distribution.

However, in radio astronomy, this assumption may no longer be valid. In order to justify the use of non-Gaussian modeling in our application, we propose in the following some numerical studies about the statistical behavior of the noise in practical scenario.

To begin with, we generate radio interferometric data according to the model in (2.21). We consider M = 8 antennas, D = 2 slighly polarized calibration sources and the associated Jones matrices are randomly simulated as Gaussian complex numbers. We study the two following cases: i) the noise n follows a typical Gaussian distribution and ii) D unknown weak sources corrupt the observations as in (2.47). We assume that all sources, calibration and non-calibration ones, are aected by similar slowly varying directional errors and we choose D = 100.

If we perform a test decision using the chi-square goodness-of-t, numerical results indicate that the data in n come indeed from a normal distribution while the realistic ones in n support the alternative hypothesis, i.e., the data do not follow such distribution at the 1% signicance level. This point is illustrated by Figures 3.1 and 3.2. We notice that the histograms for n are substantially dierent from those obtained with traditional Gaussian noise. To conrm the non-Gaussianity of the noise, we plot in Figure 3.3 the histograms of the two dierent kinds of noisy data, along with a normal density function. We also display in Figure 3.4 the quantiles of sampled data vs. the theoretical quantiles from a standard normal distribution.

We notice a linear straight line in the case of Gaussian noise n , which is consistent and deviations in the tails for n which illustrates its non-Gaussian nature. Therefore, the normal density is not so well tted when contribution of D weak unknown sources is taken into account. From the previous statistical studies, it seems clear that considering a Gaussian model for the noise in radio astronomy is not adapted.

Compound-Gaussian noise modeling

To propose a suitable alternative to the typical Gaussian noise assumption in a robust scenario, we study a wide class of distributions, gathered under the so-called CG modeling [START_REF] Jay | Détection en environnement non gaussien[END_REF][START_REF] Yao | Spherically invariant random processes: Theory and applications[END_REF].

Such distribution reveals to be suitable to achieve robustness w.r.t. outliers. Indeed, it encompasses a wide range of dierent noise distributions, among which the Student's t modeling Standard Normal Quantiles In (2.21), the noise includes contribution from all baselines, i.e., n = [n 12 , n 13 , . . . , n (M -1)M ] T .

For a particular antenna pair (p, q) ∈ {1, . . . , M } 2 , the two-scale CG distribution is generated as [START_REF] Wang | Maximum likelihood estimation of compoundgaussian clutter and target parameters[END_REF] n pq = √ τ pq g pq ,

where τ pq is a positive random variable, referred to as the texture part in the radar context, and g pq is the speckle component, given by g pq ∼ CN (0, Ω)

(3.3)
which is a zero-mean Gaussian distribution with unknown 4 × 4 speckle covariance matrix Ω.

To consolidate our choice to consider such distribution to model the noise in visibility measurements, we compare in Figure 3.5 the empirical cumulative distribution functions of the generated data n with the typical Gaussian case, the Student's t distribution as chosen in [START_REF] Kazemi | Robust radio interferometric calibration using the tdistribution[END_REF] with xed ν = 2 here and the proposed two-scale CG modeling as introduced in (3.2). We notice that the nearest curve to the realistic noise part n is obtained with the CG-based model.

Model setup

We recall that the array is made up of M antennas with known locations and the sky is composed of D calibration sources. As exposed in section 2.3, each antenna provides sensitivity to both polarization directions (x, y) and the relation between the i-th incident electromagnetic wave and the generated voltage at the p-th antenna is given in (2.9).

Non-structured case

In the so-called non-structured case, each 2 × 2 Jones matrix J i,p (θ) is parametrized by the unknown vector θ which corresponds to the entries of all Jones matrices [14,[START_REF] Nunhokee | Link between ghost artefacts, source suppression and incomplete calibration sky models[END_REF]. As we can associate a dierent Jones matrix with each signal propagation path i-p, the total number of Jones matrices is DM and thus, θ collects 8DM unknown real elements. We remind that Jones matrices model the array response and all environmental and instrumental distorsions along each propagation path.

In the radio interferometric systems under study, the measured observations consist of crosscorrelations computed for each baseline, or antenna pair (p, q) ∈ {1, . . . , M } 2 , with p < q, resulting in (2.12), if we assume uncorrelated emitted signals. The 2 × 2 source coherency matrix C i is assumed known from prior knowledge. Let us note that with the condition p < q, autocorrelations V pp are ignored which is typical in radio astronomy as they are automically agged [START_REF] Van Haarlem | LOFAR: The LOw-Frequency ARray[END_REF].

Using the vectorization operator, in the noise-free case, leads to the following 4 × 1 vector

ṽpq (θ) = vec J i,p (θ)C i J H i,q (θ) = D i=1 s i,pq (θ) (3.4)
where s i,pq (θ) = J * i,q (θ) ⊗ J i,p (θ) c i and c i = vec(C i ). For noisy measurements, we introduce a noise vector n pq for each antenna pair, leading to v pq = ṽpq (θ) + n pq (3.5) which amounts to (2.20). By stacking all possible cross-correlations in a global 4B × 1 complex vector, we obtain the full data vector y in (2.21) where s i (θ) = s T i,12 (θ), s T i,13 (θ), . . . , s T i,(M -1)M (θ) The non-structured case is particularly exible as there is no need the specify the full propagation path, thus avoiding any model misspecications, but the number of unknowns can be relatively large.

Structured case

In this section, we introduce specic models to describe the physical mechanism behind each perturbation eect: only physically meaningful parameters are estimated instead of all entries of all Jones matrices. To this end, we consider a particular context which is regime 3 of section 2.2.3 [START_REF] Lonsdale | Calibration approaches, LFD memo 015[END_REF]: sensors are clustered into a compact array and their eld-of-view are wide, resulting in direction dependent perturbations as propagation conditions dier depending on which source is targeted. But due to the small spatial extend, the same part of the ionosphere is probed and ionospheric distorsions are assumed similar for all receiving elements. Such calibration scenario is well-adapted for LOFAR or SKA stations [START_REF] Wijnholds | Calibration challenges for future radio telescopes[END_REF][START_REF] Wijnholds | Fish-eye observing with phased array radio telescopes[END_REF].

In regime 3, a particular decomposition of each Jones matrix is given by [2426, 31]

J i,p (θ i,p ) = G p (g p )H i,p Z i (ϕ i )F i (ϑ i ) (3.6)
for i ∈ {1, . . . , D}, p ∈ {1, . . . , M } and θ i,p = [ϑ i , ϕ i , g T p ] T . • Ionospheric eects: As described in section 2.2.1, the ionosphere is a disturbed environment with spatial variations. One of the prevailing eects is a propagation delay which results in apparent shift of the source position [START_REF] Wijnholds | Calibration challenges for future radio telescopes[END_REF][START_REF] Van Der Tol | Self-calibration for the LOFAR radio astronomical array[END_REF]. We adapt (2.16) to our specic scenario, leading to

Z i (ϕ i ) = exp jϕ i I 2 . (3.7)

CHAPTER 3. ROBUST CALIBRATION METHOD IN NON-GAUSSIAN ENVIRONMENT

Propagation through the ionosphere also results in rotation eects whose form is given in (2.18). For the ionospheric Faraday rotation, we write it as

F i (ϑ i ) = cos(ϑ i ) -sin(ϑ i ) sin(ϑ i ) cos(ϑ i ) (3.8)
where ϑ i is the unknown Faraday rotation angle.

As the array is compact with similar elements in regime 3, propagation delay (3.7) and rotation (3.8) are assumed identical for all antennas here [START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF].

• Instrumental eects: Each antenna is described by two complex electronic gains, s.t.,

G p (g p ) = diag{g p } (3.9)
where g p is the 2 × 1 unknown complex gain vector.

Finally, H i,p is an assumed known matrix for which prior information is provided by the available calibration source locations and antenna positions in (2.15) but also thanks to electromagnetic simulations as regards the antenna response and the beam pattern [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF][START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF][START_REF] Yatawatta | Reduced ambiguity calibration for LOFAR[END_REF][START_REF] Kazemi | Blind calibration for radio interferometry using convex optimization[END_REF].

In this structured case, unknown parameters of interest are given by θ i,p for i ∈ {1, . . . , D} and p ∈ {1, . . . , M }. We collect them all in ε = P[θ T 1,1 , θ T 1,2 , . . . , θ T D,M ] T where P is an appropriate rearrangement matrix, avoiding redundancy and leading to

ε = [ϑ 1 , . . . , ϑ D , ϕ 1 , . . . , ϕ D , g T 1 , . . . , g T M ] T . (3.10)
Therefore, the global number of real unknowns is 2(D + M ) and we assume M D so that observations outnumber unknown parameters.

Estimation of the texture and speckle parameters

Estimation of the noise parameters in (3.2) can be performed in two ways: either a priori knowledge is available about the distribution of τ pq and a Bayesian approach is conducted or the distribution is unknown and it can be assumed deterministic during the estimation procedure. All the algorithms we are proposing in this work are iterative and based on the ML estimator [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF][START_REF] Ollier | Joint ML calibration and DOA estimation with separated arrays[END_REF]

Bayesian approach

Iterative MAP estimator (IMAPE): From a Bayesian perspective, we investigate robust calibration considering assumed specic distributions for the texture realizations τ pq , resulting in MAP-based estimators [START_REF] Ollier | Bayesian calibration using dierent prior distributions: an iterative maximum a posteriori approach for radio interferometers[END_REF]. Using Bayes' theorem, the expression of the joint log-likelihood function, denoted as L J , is given by

L J = L C + pq ln p(τ pq ; ϕ) (3.11)
where L C is the conditional likelihood function and ϕ refers to the unknown hyperparameters describing the corresponding prior distribution. Thus, such parameters need to be estimated as well. The principle of the MAP estimator is to maximize (3.11) w.r.t. each unknown individual parameter through a step-wise approach, leading to alternative update of unknowns θ, ϕ, Ω and τ . Assuming independence between n pq , the expression of L C is given by As it can been seen from (3.2), the probability density function (pdf ) of each texture parameter τ pq needs to be specied. Depending on the nature of τ pq , dierent algorithms can be developed. We exploit prior distributions for the texture, available in closed-form, which generate various heavy-tailed noise models for n pq . In the following, we dene specic texture priors and derive the corresponding IMAPE for each prior distribution p(τ pq ; ϕ).

L C = ln p (y|τ ; θ, Ω) = -
• Gamma distribution: The corresponding pdf is given by

p(τ pq ; a, b) = 1 Γ(a)b a τ a-1 pq exp - τ pq b (3.13)
where a and b are the shape and scale parameters, respectively. Such prior leads to the Kdistribution for the noise vector n pq [START_REF] Ollila | Complex elliptically symmetric distributions survey, new results and applications[END_REF].

If we insert (3.13) into (3.11), the resulting joint log-likelihood function reads

L J = L C + (a -1) pq ln(τ pq ) - pq τ pq b -B ln Γ(a) -Ba ln(b).
Considering ∂L J /τ pq = 0 leads to the following estimate for the texture parameter which can be computed numerically.

τ pq = (a -5)b + (a -5) 2 b 2 + 4bu H pq (θ)Ω -1 u pq (θ)
In this specic case, steps of the corresponding IMAPE are described in Algorithm 1. Due to the iterative approach, we notice that (3.15) was plugged into (3.16) as well as (3.14) into (3.28), see section 3.3.2. Similar algorithms can be derived for each texture prior.

• Inverse Gamma distribution: When the texture parameter follows an inverse Gamma function, s.t.,

p(τ pq ; a, b) = b a Γ(a) τ -a-1 pq exp - b τ pq , (3.17)
the generated noise model is a Student's t distribution [START_REF] Lange | Robust statistical modeling using the t distribution[END_REF][START_REF] Ollier | Bayesian lower bounds for dense or sparse (oulier) noise in the RMT framework[END_REF]. Let us note that when the shape parameter is xed in (3.17), s.t., a = 1, the CG-based model for n pq becomes a Cauchy distribution [START_REF] Jay | Détection en environnement non gaussien[END_REF].

Using (3.17), (3.11) and the derivative w.r.t. τ pq , the corresponding estimate reads We mentioned the particular case of Cauchy distribution for which (3.20) becomes â = 1.

τ pq = b + u H pq (θ)Ω -1 u pq (θ) a + 5 .
• Exponential distribution: Such prior for the texture reads p(τ pq ; λ) = λ exp (-λτ pq ) , (3.21) where λ is the so-called rate parameter. The resulting model for the noise is a Laplace distribution [START_REF] Jay | Détection en environnement non gaussien[END_REF].

Based on (3.21) and (3.11), the calculations result in

τ pq = -4 + 16 + 4λu H pq (θ)Ω -1 u pq (θ) 1/2 2λ (3.22)
and the single rate parameter is deducted from λ = B pq τ pq .

(3.23)

• Inverse Gaussian distribution: The inverse Gaussian distribution, also called Wald distribution, with shape λ and assumed unit mean [START_REF] Ollila | Complex elliptically symmetric distributions survey, new results and applications[END_REF], is described by the following pdf:

p(τ pq ; λ) = λ 2π 1/2 τ -3/2 pq exp -λ(τ pq -1) 2 2τ pq . (3.24)
The resulting CG model is referred to as the IG-CG distribution.

Estimates of noise parameters are given by

τ pq = -11 + 121 + 4λ λ + 2u H pq (θ)Ω -1 u pq (θ)
To sum up, we expose in Table 3.1 the dierent kinds of CG-based noise models we are considering in this study and their corresponding texture priors with estimates.

τ pq φ Gamma p(τ pq ; a, b) = 1 Γ(a)b a τ a-1 pq exp - τ pq b K (3.14) (3.16) & (3.15) Inverse Gamma p(τ pq ; a, b) = b a Γ(a) τ -a-1 pq exp -b τ pq , Student's t (3.18) (3.20) & (3.19) Inverse Gamma p(τ pq ; b) = bτ -2 pq exp -b τ pq , Cauchy (3.18)
a = 1 & (3.19) Exponential p(τ pq ; λ) = λ exp (-λτ pq ) Laplace (3.22) (3.23) Inverse Gaussian p(τ pq ; λ) = λ 2π 1/2 τ -3/2 pq exp -λ(τ pq -1) 2 2τ pq IG-CG (3.25) (3.26)

(Relaxed) deterministic approach

The noise model is not exactly known in practice, so prior knowledge about the pdf is not always available. Thus, the statistical distribution of the texture parameter being unknown, we can consider it as an unknown deterministic parameter in the estimation process [START_REF] Conte | Recursive estimation of the covariance matrix of a compound-gaussian process and its application to adaptive CFAR detection[END_REF][START_REF] Pascal | Covariance structure maximum-likelihood estimates in compound gaussian noise: Existence and algorithm analysis[END_REF]. This avoids to choose a specic distribution which may lead to model misspecications and enables to propose a broad robust estimator w.r.t. the presence of outliers. We call it the relaxed deterministic approach as it is not strictly based on the exact model in (3.2) [START_REF] Ollier | Relaxed concentrated MLE for robust calibration of radio interferometers[END_REF][START_REF]Robust calibration of radio interferometers in non-Gaussian environment[END_REF].

Relaxed concentrated ML estimator (RCMLE): Optimization of (3.12) is performed w.r.t.

each unknown parameter sequentially, while the others are xed. Closed-form expressions can be obtained for the texture component τ pq , for which no distribution has to be specied, and the speckle covariance matrix Ω. Estimation of θ is the step of interest and will be discussed in section 3.4. If we take the derivative of (3.12) w.r.t. τ pq and equate it to 0, we obtain the following expression for the texture estimate

τ pq = 1 4 u H pq (θ)Ω -1 u pq (θ). (3.27)
Using classical dierential properties [140, p. 2741] and the permutation property, estimation of the speckle covariance matrix reads

Ω = 1 B pq 1 τ pq u pq (θ)u H pq (θ). (3.28)
We adopt an iterative procedure here with a concentrated ML scheme. Thus, by plugging (3.27) into (3.28), we nally obtain

Ωt+1 = 4 B pq u pq (θ)u H pq (θ) u H pq (θ) Ωt -1 u pq (θ) (3.29)
with t the iteration counter. Let us note that a constraint is required to remove scaling ambiguities in model (3.2). Therefore, we impose tr {Ω} = 1 but this choice is arbitrary and does not aect the estimates of interest [START_REF] Zhang | MIMO radar target localization and performance evaluation under SIRP clutter[END_REF]. The normalization requires the following step in the procedure Ωt+1 = Ωt+1 tr Ωt+1 .

(3.30)

The global scheme of the proposed RCMLE is described in Algorithm 2.

Algorithm 2: RCMLE Relaxed concentrated ML estimator input : D, M , B, {C i } i∈{1,...,D} , y output : Remark: Let us note that it is possible to consider a dierent covariance matrix Ω pq for each baseline (p, q) in (3.3). In this specic case, the proposed robust calibration algorithm requires a few modications and corresponding expressions are presented in Appendix A.

θ initialize: Ω ← Ω init , τ ← τ init , θ ←

Estimation of Jones matrices

The proposed robust calibration estimator is based on the CG model in (3.2). We assume that noise parameters τ pq and Ω are estimated with the Bayesian or deterministic approach.

For sake of clarity, we consider the relaxed concentrated ML method for which texture realizations are assumed deterministic. But the corresponding expressions for the Bayesian case are straightforward to obtain, using section 3.3.1.

Specically, each block of unknown parameters is updated sequentially in a global iterative loop as shown in Algorithm 2. Estimates of noise parameters, i.e., the texture and the speckle parts, are derived from (3.27) and (3.29), respectively, followed by the normalization constraint in (3.30), for identiability issues.

Non-structured case

The goal is to estimate θ, considering (3.12). For a given Ω and τ , it consists in the following minimization problem θ = argmin

θ pq 1 τ pq u H pq (θ)Ω -1 u pq (θ) . (3.31)
Thus, a multi-dimensional optimization approach needs to be carried out to solve (3.31), resulting in signicant computation time due to the large size of θ. To reduce the computational cost, we propose to use the EM algorithm [START_REF] Moon | The Expectation-Maximization algorithm[END_REF][START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF] which requires a proper parameter partitioning and the denition of complete data. In the non-structured case, it is natural to decompose the parameter of interest θ as follows

θ = [θ T 1 , . . . , θ T D ] T = [θ T 1,1 , . . . , θ T 1,M , . . . , θ T D,1 , . . . , θ T D,M ] T .
(3.32)

ESTIMATION OF JONES MATRICES 43

This way, each propagation path i-p is parametrized by the vector θ i,p ∈ R 8×1 , s.t., J i,p (θ) = J i,p (θ i,p ).

As explained in section 2.4.3.1, the EM algorithm iterates between two steps and provides the ML estimates with limited complexity. These two steps are explained in the following

• E-step: The so-called complete data vector is given by (2.44) with a decomposition into source contributions s.t. y = D i=1 w i and n = D i=1 n i . Each noise part is generated as

n i ∼ CN (0, β i Ψ) (3.33)
Since measurements are assumed independent between antenna pairs, as seen in (3.12), the covariance Ψ has the following block-diagonal expression Ψ = bdiag{τ 12 Ω, . . . , τ (M -1)M Ω}. (3.37)

• M-step: Instead of directly estimating θ ∈ C 4DM ×1 from (3.31), which can be computationally prohibitive, the global multiple source estimation problem is reduced to multiple single source sub-problems. Indeed, once ŵ is evaluated by the previous E-step, we need to maximize the following likelihood function

p( ŵ|θ, τ , Ω) = 1 |πΞ| exp -ŵ -s(θ) H Ξ -1 ŵ -s(θ) = D i=1 1 |πβ i Ψ| exp -ŵi -s i (θ i ) H (β i Ψ) -1 ŵi -s i (θ i ) (3.38)
with assumed independent w i . Therefore, estimation of θ i ∈ C 4M ×1 is obtained through minimization of the following objective function

φ i (θ i ) = ŵi -s i (θ i ) H (β i Ψ) -1 ŵi -s i (θ i ) . (3.39)
This can be performed numerically with the LM algorithm [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Madsen | Methods for non-linear least squares problems[END_REF][START_REF] Gavin | The Levenberg-Marquardt method for nonlinear least squares curve-tting problems[END_REF], as exposed in section 2.4.2.1. However, obtaining closed-form expressions enables to reduce even more the complexity of the problem. To this end, we apply the BCD algorithm [START_REF] Friedman | Pathwise coordinate optimization[END_REF][START_REF] Hong | A unied algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing[END_REF] for each single source subproblems in (3.39).

As shown in (3.32), the parameter vector is partitioned according to each source but also to each antenna. Thus, in the following, instead of optimizing (3.39) w.r.t. θ i , we propose to carry out the optimization w.r.t. θ i,p , for xed θ i,q with q = p. We will see that an analytical solution is obtained for θi,p as function of all other θ i,q with q = p. The procedure is repeated alternatively for each component vector θ i,p with p ∈ {1, . . . , M } until convergence. In (3.39), only s i is function of θ i but actually, only a subset depends on the block-coordinate vector θ i,p , i.e., {s i,pq } for q > p, q ∈ {1, . . . , M } and {s i,qp } for q < p, q ∈ {1, . . . , M }. Thus, it is more convenient to rewrite (3.39) as follows

φ i (θ i,p ) = M q=1 q>p w i,pq -s i,pq (θ i,p ) H (β i τ pq Ω) -1 w i,pq -s i,pq (θ i,p ) + M q=1 q<p w i,qp -s i,qp (θ i,p ) H (β i τ qp Ω) -1 w i,qp -s i,qp (θ i,p ) + Constant. (3.40)
In the Constant part, we incorporate all subsets of φ i (θ i,p ) which do not depend on θ i,p . Notation, calculus and details are given in Appendix B in which we show that it is possible to write

s i,pq (θ i,p ) = Σ i,q θ i,p (3.41)
and, likewise

s i,qp (θ i,p ) = Υ i,q θ * i,p . (3.42)
In the end, the expression of θ i,p for p ∈ {1, . . . , M } and the given i-th source direction is written

as θi,p =    (Σ H i A i,p Σ i + Υ H i Ãi,p Υ i ) -1 (Σ H i A i,p w i,p + Υ H i Ãi,p wi,p ) for 1 < p < M (Σ H i A i,p Σ i ) -1 Σ H i A i,p w i,p for p = 1 (Υ H i Ãi,p Υ i ) -1 Υ H i Ãi,p wi,p for p = M (3.43)
Finally, from the above discussion, the detailed scheme of the proposed RCMLE is summarized in Algorithm 3.

Algorithm 3: RCMLE Relaxed concentrated ML estimator input : D, M , B, {C i , β i } i∈{1,...,D} , y output :

θ initialize: Ω ← Ω init , τ ← τ init , θ ← θ init
while stop criterion unreached (referred to as the rst loop) do while stop criterion unreached (referred to as the second loop) do 1 E-step: ŵi obtained from (3.37), i ∈ {1, . . . , D} Regarding the convergence properties of the proposed algorithm, we notice that: the maximization step of the EM algorithm is computed thanks to the BCD algorithm by considering the cost function in (3.39). Minimization leads to a unique solution for each individual block of variables, given in (3.43). Thus, convergence to a stationary point and local convergence are ensured [START_REF] Bertsekas | Nonlinear programming[END_REF]. When the M-step is solved exactly, i.e., when the BCD algorithm provides the exact minimizer of (3.39), and for a theoretical innite number of iterations, convergence to a stationary point by the EM algorithm is ensured. Still, let us note that depending on the initialization, the stationary point can be either a global, a local minimizer or, in some unusual cases, a saddle point [START_REF] Mclachlan | The EM algorithm and extensions[END_REF]. Proper initialization is thus required based on possible prior information or one can use some metaheuristic approaches such as random-restart hill climbing or SA methods. If the EM provides exact solutions, monotone local convergence of the ML estimator is ensured since the value of the cost function at each step can either improve or maintain but cannot worsen [START_REF] Vorobyov | Maximum likelihood direction-of-arrival estimation in unknown noise elds using sparse sensor arrays[END_REF]. Therefore, convergence to a local stationary point is ensured for the overall algorithm if convergence is attained in each loop (three loops as shown in Algorithm 3), which requires an adequate initialization and a theoretical innite number of iterations. However, in practice, numerical simulations in section 3.5 show that a few iterations are enough to reach convergence to, at least, a local minimizer and achieve relatively good numerical stability.

Structured case

In the structured case, we consider the specic model in (3.6), adapted to regime 3. The RCMLE in Algorithm 3 provides an estimate for all Ĵi,p with i ∈ {1, . . . , D} and p ∈ {1, . . . , M }.

The idea is to estimate the unknown parameter vector of interest ε in a sequential manner, once RCMLE has been performed and all Jones entries have been estimated. To this end, we adopt an iterative procedure and optimize a given cost function w.r.t. one of the physical parameters in ε while xing the others, leading to a global alternative method.

• Estimation of g p : Electronic gains are obtained by solving ĝp = argmin gp κ(g p )

(3.44) in which κ(g p ) = D i=1 || Ĵi,p -G p (g p )H i,p Z i F i || 2 F = D i=1 Tr Ĵi,p -G p (g p )R i,p Ĵi,p -G p (g p )R i,p H (3.45)
where R i,p = H i,p Z i F i . Taking the derivative w.r.t.

[g p ] k for k ∈ {1, 2} leads to ∂κ(g p ) ∂[g p ] k = D i=1 Tr -e k e T k R i,p ĴH i,p + e k e T k R i,p R H i,p G H p . (3.46)
Equating (3.46) to 0 results in D i=1 [X i,p ] k,k = D i=1 [W i,p ĜH p ] k,k = D i=1 [W i,p ] k,k [ Ĝ * p ] k,k (3.47)
where X i,p = R i,p ĴH i,p and W i,p = R i,p R H i,p . Each sensor gain element is thus estimated by

[ĝ p ] k = D i=1 [W * i,p ] k,k -1 D i=1 [X * i,p ] k,k . (3.48)
• Estimation of ϕ i : To do so, we consider the following minimization problem φi = argmin

ϕ i κ(ϕ i ) (3.49) where κ(ϕ i ) = M p=1 || Ĵi,p -G p H i,p Z i (ϕ i )F i || 2 F .
Considering ∂κ(ϕ i )/∂ϕ i and setting the result to zero leads to

M p=1

Tr j exp -j φi Ĵi,p

F H i H H i,p G H p -j exp j φi G p H i,p F i ĴH i,p = 0 (3.50)
and nally

exp 2j φi = Tr M i,p Tr M H i,p (3.51)
with M i,p = M p=1 Ĵi,p F H i H H i,p G H p and from which we can directly deduce φi . • Estimation of ϑ i : We address the following one-dimensional minimization problem θi = argmin

ϑ i M p=1 || Ĵi,p -G p H i,p Z i F i (ϑ i )|| 2 F (3.52)
for each source. Estimates can be obtained in a reasonable computational time through classical data grid search followed by Newton type algorithm.

The global procedure of the RCMLE in the structured case is exposed in Algorithm 4.

Simulation results

In this section, we aim to evaluate the performance of the designed algorithm: we start with a CG-based scenario where the noise model matches our noise assumption and we plot the corresponding CRB. Afterwards, we intend to compare, in a more realistic scenario, the proposed algorithm with the recently introduced approach based on the Student's t [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF] and the traditional Gaussian cases [14]. We study both non-structured and structured situations and also compare the relaxed deterministic technique with Bayesian alternatives. Finally, visual performance is provided on academic examples by creating residual images where weak background sources are revealed.

Under CG-based noise

Firstly, we aim to assess the estimation performance of the proposed CG-based RCMLE and its numerical stability. We consider the non-structured case where θ in (3.32) corresponds to the real and imaginary parts of the entries of all Jones matrices. We assume that the additive noise in (3.5) follows a CG distribution as in (3.2), i.e., the generated noise matches the noise model assumption. We choose to generate each random texture component from an inverse Gamma distribution [START_REF] Balleri | Maximum likelihood estimation for compound-Gaussian clutter with inverse gamma texture[END_REF], as follows

τ pq ∼ IG(ν/2, ν/2), (3.53)
with ν degrees of freedom [START_REF] Ollila | Complex elliptically symmetric distributions survey, new results and applications[END_REF] and [Ω] k,l = σ 2 0.9 |k-l| exp j π 2 (k-l) . With (3.53), the resulting noise distribution for n pq is a Student's t.

To assess the statistical performance, we make use of the CRB [START_REF] Stoica | Spectral analysis of signals[END_REF]. Under quite general/weak conditions, the variance satises

MSE([ θ] k ) = E [ θ] k -[θ] k 2 [CRB(θ)] k,k (3.54)
where the CRB is the inverse of the FIM, noted F. In our specic case, we adapt the Slepianbangs type formula to our context [START_REF] Besson | On the Fisher information matrix for multivariate elliptically contoured distributions[END_REF][START_REF] El Korso | CRLB under K-distributed observation with parameterized mean[END_REF] and obtain

[F] k,l = 2 ν + 4 ν + 5 pq ∂ ṽH pq (θ) ∂[θ] k Ω -1 ∂ ṽpq (θ) ∂[θ] l . (3.55)
Noise parameters being decoupled from parameters of interest θ, we only keep the part of the FIM corresponding to the latter.

We consider D = 2 calibration sources and M = 8 sensors in the array, so the total number of real unknowns is given by 8DM = 128 and the number of measurements is 8B = 224. Prior information about the source C i is generated thanks to random Stokes parameters [START_REF] Hamaker | Understanding radio polarimetry[END_REF][START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF]. In

MSE = 1 M C M C k=1 ( θk -θ) 2 (3.56)
where M C is the number of Monte-Carlo runs, which we usually x to 100 in our work. In Indeed, the RCMLE derives estimates for unknown and deterministic texture parameters while they were in fact generated as inverse Gamma distributed random parameters (3.53) in the data model and the CRB relies on this specic prior distribution.

In Figure 3.7, we investigate numerical convergence rates of the RCMLE, for each of its three loops as exposed in Algorithm 3, and see how estimation of θ is aected through time. To do so, we introduce the following quantity

t {θ} = || θ t -θ t-1 || 2 2 .
(3.57)

Under realistic scenario

We aim to study the designed algorithm in a realistic scenario with D calibration sources and D weak background sources which act as outliers and contribute to the noise component, as shown in (2.47) and (3.1). The SNR is dened as the ratio of the normalized power of D calibration sources over the sum of normalized power of D background sources and a power noise factor. Jones matrices stand for the perturbation eects along the signal propogation path and robustness. Besides, no assumption has been made about independent entries within the noise vector n pq , unlike [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF]. Let us remark that in Figure 3.9(a), the MSE is plotted for both D = 2 calibration sources, leading to 4DM = 64 plotted points for the real part of θ.

Under structured case

Once all Jones entries in θ have been estimated, we aim to apply the algorithm for structured Jones matrices, exposed in Algorithm 4. Both algorithms based on the Student's t [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF] and the Gaussian distribution [14] have been introduced in the non-structured case so we apply the structured case on the output of these two approaches. Jones matrices are randomly generated according to (3.6) and we introduce g = [g T 1 , . . . , g T M ] T . For the ionospheric phase delay in (3.7), it is generated as follows [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF] ϕ i,p = η i u p + ζ i v p .

(3.58)

3.5. SIMULATION RESULTS

51

Equivalently, we can write

ϕ T i = α T i Λ (3.59)
where

ϕ i = [ϕ i,1 , . . . , ϕ i,M ] T , Λ = u 1 , . . . , u M v 1 , . . . , v M , α i = [η i , ζ i ]
T refers to the two-dimensional osets by which the i-th source is shifted and known position r p = [u p , v p ] T of the p-th antenna is expressed in wavelength units. Once ϕ i,p has been estimated with (3.51) for p ∈ {1, . . . , M }, we can deduce directional shifts αi with the following aditional step

αT i = φT i Λ H M p=1 v 2 p -M p=1 u p v p -M p=1 v p u p M p=1 u 2 p M p=1 u 2 p M p=1 v 2 p -( M p=1 u p v p ) 2
.

(3.60)

The MSE results are shown in Figure 3.10, with similar computation times, for specic physical parameters: the complex gains for xed SNR and the source oset ζ 1 as function of the SNR, the behavior being the same for other parameters. Better performances are obtained with the proposed RCMLE in the structured case. This was already expected from Figure 3.9 as better estimation of Jones entries naturally leads to better estimation of parameters describing them.

Recovered images

In this section, we investigate the formation of residual images on academic examples, obtained with Matlab, but more realistic data simulation and image recovering will be presented in Chapter 5.

Let us consider D = 2 calibration sources in the sky model. In the ideal case, visibilities are only aected by the K-Jones term in (2.15). But due to perturbation eects, especially the ionospheric phase delays in (3.58), source positions are subject to shifts and the apparent location may dier from the exact one. We wish to recover D = 4 background sources whose exact intensities and positions are unknown to us and randomly generated from a discrete uniform distribution. Only ionospheric phase delays are incorporated within the Jones chain (2.13).

Calibration is performed with three dierent methods: the proposed RCMLE, the Student's t [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF] and the Gaussian cases [14]. From the phase estimates, we deduce the osets with (3.60) by which the D = 2 calibration sources are shifted. Then, the corresponding osets for weak sources are deduced by linear interpolation. More specically, we are interested in residual visibilities. Indeed, from (2.47), for each antenna pair, we know that

v pq = D i=1 s i,pq (θ) + D i =1 s i ,pq + n pq (3.61)
with the calibration, the non-calibration and the noise components, respectively. By residual visibilities, we mean the following quantity

v RES pq = v pq - D i=1 s i,pq ( θ) (3.62)
in which θ is estimated using one of the three aformentionned estimators. The true situation corresponds to the case when θ is exactly known. Intensities of weak sources are then recovered by performing a two-dimensional Inverse Fourier Transform (2.6) on the residual visibilities. 3.11 for a mostly empty sky with point-like sources. Thus, if calibration is performed exactly, there should be no ux left at the calibration positions. Position shifts being hardly visible, to make it more clear, we compute in Table 3.2 the mean over all D sources of the dierence between the exact source location and the estimated one, in both directions (x, y).

The lower is this dierence, the more accurate is the calibration method. We recover similar residual images in Figure 3.12 but intensity is shown for each position of the plane.

Conclusion

In radio astronomy, Jones matrices model all the perturbation eects which corrupt the observations and they can be considered as structured or not depending on the scenario under study. However, visibility measurements are also aected by the presence of outliers wich change the typical Gaussian nature of noise. To achieve robustness, we propose the RCMLE based on the CG model and a relaxed assumption. This algorithm reveals to be more robust for both non-structured and structured cases, it is less subject to errors than MAP approaches and benets from reasonable computational complexity. Indeed, using the EM and the BCD algorithms enables to reduce the global computational cost thanks to a per source and per antenna parameter partition, leading to multiple sub-problems instead of a global multi-dimensional one.

Texture, speckle components and Jones matrices are estimated iteratively thanks to closed-form expressions. However, their specic structure of variation, notably w.r.t. frequency, is not taken into account and could be exploited to enhance the calibration process. Thus, the extension to the multi-frequency scenario is studied in the following chapter.

Variation of parameters across frequency

We are now dealing with interferometer calibration where direction independent and direction dependent eects need to be corrected. The Jones chain given in (3.6) for regime 3 is adapted to regime 4 in a multi-frequency scenario as follows

J [f ] i,p (θ [f ] ) = G [f ] p (g [f ] p )H [f ] i,p Z [f ] i,p (ϕ [f ] i,p )F [f ] i,p (ϑ [f ] i,p ) (4.1)
in which the frequency dependence of ionospheric distortions can be accurately described. Indeed, the ionospheric phase delay in (3.7) is specied now as [START_REF] Van Der Tol | Ionospheric calibration for the LOFAR radio telescope[END_REF] ϕ

[f ] i,p ∝ T EC i,p f (4.2)
where T EC was dened in section 2.2.1. As regards the Faraday rotation in (3.8), it is given by [START_REF] Smirnov | Revisiting the radio interferometer measurement equation[END_REF] ϑ

[f ] i,p ∝ RM i,p f 2 (4.3)
in which RM refers to the rotation measure, depending on the magnetic eld and the density of electrons along the propagation path i-p.

There is no available model for the frequency variation of the gains g

[f ] p . However, we can enforce smoothness across frequency by introducing polynomial variation across frequency, as discussed in section 4.3.2.

Therefore, for one given frequency f ∈ F = {f 1 , . . . , f F }, the number of unknown parameters of interest is gathered into the (2M

D + 2M ) × 1 vector θ [f ] = [ [f ] T , g [f ] T 1 , . . . , g [f ] T M ] T where [f ] = [ϑ [f ] 1,1 , . . . , ϑ [f ] D,M , exp(jϕ [f ]
1,1 ), . . . , exp(jϕ (M -1)M ] T and the 4 × 4 speckle covariance matrix Ω [f ] . We notice that frequency dependence has been added for the noise parameters as well but we do not assume any specic prior structure w.r.t. f . Nevertheless, if a priori knowledge about noise variation w.r.t. frequency is available, it can be straightforwardly incorporated into the proposed algorithm.

4.3

Estimation of Jones matrices in structured case

Robust estimation in multi-frequency scenario

Assuming independence between n

[f ] pq , the corresponding expression of the log-likelihood in (3.12) reads

ln p {y [f ] } f ∈F | {θ [f ] , τ [f ] , Ω [f ] } f ∈F = - f ∈F l [f ] (θ [f ] ) - f ∈F pq ln |πτ [f ] pq Ω [f ] | (4.4) with l [f ] (θ [f ] ) = pq 1 τ [f ] pq u [f ] H pq (θ [f ] )Ω [f ] -1 u [f ] pq (θ [f ]
). Then, we adopt the same procedure as for the RCMLE and obtain similar texture and speckle covariance estimates than in (3.27) and (3.29). More specically, we can write

τ [f ] pq = 1 4 u [f ] H pq (θ [f ] )Ω [f ] -1 u [f ] pq (θ [f ] ) (4.
Ω[f] = 4 B pq u [f ] pq (θ [f ] )u [f ] H pq (θ [f ] ) u [f ] H pq (θ [f ] )( Ω[f] ) -1 u [f ] pq (θ [f ] ) (4.7)
for the speckle, followed by the normalization step

Ω[f] ← Ω[f] tr Ω[f]
to avoid any ambiguity.

The principle of the proposed multi-frequency robust approach is similar to the one exposed in Algorithm 2: we estimate alternatively { θ

[f] } f ∈F , { Ω[f] } f ∈F and {τ [f ] } f ∈F in a global iterative loop.
In what follows, we focus on estimation of {θ [f ] } f ∈F , for a xed {τ [f ] } f ∈F given by (4.6) and { Ω[f] } f ∈F given by (4.7), by considering consensus-based distributed optimization [START_REF] Bertsekas | Parallel and distributed computation: numerical methods[END_REF][START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF].

Distributed estimation in multi-frequency scenario

Multi-frequency calibration considers a whole frequency range with multiple sub-frequency bands. A computationally ecient way to handle this situation is to apply decentralized strategies with distributed and consensus algorithms as the ADMM which was introduced in section 2.4.3.4. This technique is well-suited for large-scale problems and has already been applied for image reconstruction [START_REF] Ferrari | Distributed image reconstruction for very large arrays in radio astronomy[END_REF]. We will be considering this procedure in the following with a group of computational agents. We assume that each of them has access to data for a specic frequency f ∈ F = {f 1 , . . . , f F } and solves a local problem. Estimates are then transmitted to the fusion center which enforces consensus among agents and nally transfers the updates to each local agent. Compared to independent (per-channel) calibration performed for one frequency after another, distributed calibration can lead to better accuracy since information across all frequencies is taken into account.

The idea is to be able to get closed-form expressions for the parameters of interest. To this end, for i ∈ {1, . . . , D} and p ∈ {1, . . . , M }, we rewrite the ionospheric phase delay as exp jϕ

[f ] i,p = b [f ] T z i,p (4.8)
and the Faraday rotation angle as

ϑ [f ] i,p = 1 f 2 z i,p (4.9)
s.t. the (M + N M ) × 1 unknown frequency independent vector of hidden variables z i = [z i,1 , . . . , z i,M , z T i,1 , . . . , z T i,M ] T is introduced for the i-th source direction. The frequency dependent vector b

[f ] is dened s.t. [b [f ] ] k = 1 f k-1 (4.10)
for k ∈ {1, . . . , N } and (N -1) is the order of the truncated approximation power series of exp jϕ

[f ] i,p , chosen arbitrarily. More specically, we can write exp jϕ

[f ] i,p N k=1 1 f k-1 (jϕ [f 0 ] i,p ) k-1 (k -1)! (4.11) s.t. ϕ [f ] i,p = ϕ [f 0 ] i,p f and ϕ [f 0 ] i,p ∝ T EC i,p . Let us note that z i,p ∝ RM i,p .
As mentioned before, we impose smoothness over frequency for the gains, s.t.,

g [f ] p = B[f] zp (4.12)
where the 2 × 2 Ñ frequency modeling matrix reads

B[f] = b[f] T ⊗ I 2 (4.13)
and the Ñ ×

1 vector b[f] is described by [ b[f] ] k = f -f 0 f 0 k-1 (4.14)
with k ∈ {1, . . . , Ñ }, an arbitrarily chosen Ñ and a given reference frequency f 0 [START_REF] Yatawatta | Distributed radio interferometric calibration[END_REF]. Frequency independent variables are introduced as well for the p-th receiving element through the 2 Ñ × 1 vector zp .

With (4.4), (4.8), (4.9) and (4.12), the calibration problem can be formulated as

{ θ[f] } f ∈F , ẑ = argmin {θ [f ] } f ∈F ,z f ∈F l [f ] (θ [f ] ) s.t. θ [f ] = B [f ] z (4.15) with z = [z 1,1 , . . . , z D,M , z T 1,1 . . . , z T D,M , zT 1 , . . . , zT M ] T the (DM + DN M + 2 Ñ M ) × 1 augmented
vector of hidden variables and the frequency coupling matrix is given by

B [f ] = bdiag 1 f 2 I DM , I DM ⊗ b [f ] T , I M ⊗ B[f] =    1 f 2 I DM 0 DM ×DN M 0 DM ×2 Ñ M 0 DM ×DM I DM ⊗ b [f ] T 0 DM ×2 Ñ M 0 2M ×DM 0 2M ×DN M I M ⊗ B[f]    . (4.16)
Therefore, multi-frequency calibration amounts to solve a global constrained optimization problem. To achieve this goal, we introduce the following Lagrangian

L {θ [f ] } f ∈F , z, {x [f ] } f ∈F = f ∈F L [f ] θ [f ] , z, x [f ] (4.17) with L [f ] θ [f ] , z, x [f ] = l [f ] θ [f ] + h [f ] θ [f ]
, z, x [f ] where the penalty term reads

h [f ] θ [f ] , z, x [f ] = 2 x [f ] H θ [f ] -B [f ] z + ρ||θ [f ] -B [f ] z|| 2 2 (4.18)
The associated Lagrange parameters (or dual variables) for a given frequency f are denoted as the (2DM + 2M )

× 1 vector x [f ] = [x [f ] 1,1 , . . . , x [f ] D,M , x [f ] 1,1 , . . . , x [f ] D,M , x[f] T 1 , . . . , x[f] T M ]
T and ρ is a regularization factor.

To sum it up, the ADMM algorithm consists in updating the three following quantities

• θ[f] t+1 = argmin θ [f ] L [f ] θ [f ] , (ẑ) t ,
• (ẑ) t+1 = argmin z f ∈F h [f ] θ[f] t+1 , z, x[f] t
performed globally at the fusion center (4.20)

• x[f] t+1 = x[f] t + ρ θ[f] t+1 -B [f ] (ẑ)
ẑ =   ρ f ∈F B [f ] H B [f ]   -1   f ∈F B [f ] T x [f ] + ρθ [f ]   (4.22)
which is then broadcasted to each agent as a common global variable in order to perform the update in (4.21) locally.

The step of interest is (4.19) and we will be using an iterative approach to estimate θ [f ] .

On the one hand, we consider [f ] and minimize successively w.r.t. each ϑ [f ] i,p and ϕ

[f ] i,p for i ∈ {1, . . . , D} and p ∈ {1, . . . , M }. On the other hand, we will be considering g [f ] . Specically,

• Faraday rotation: Estimation of ϑ [f ] i,p amounts to solve θ[f] i,p = argmin ϑ [f ] i,p L [f ] θ [f ] , z, x [f ] (4.23)
which can be computed with a one-dimensional Newton or gradient descent-type algorithm [START_REF] Nocedal | Numerical optimization[END_REF].

Parallelization is also possible for all i ∈ {1, . . . , D} and p ∈ {1, . . . , M }. i,p results, after some calculus, in

φ[f] i,p = 1 2 arg - α [f ] i,p β [f ] i,p . (4.24)
Notations and detailed calculations are provided in Appendix C.

• Electronic gains:

In what follows, we omit dependence w.r.t. parameters of interest for sake of clarity. Let us recall that g

[f ] p = [g [f ] p ] 1 , [g [f ] p ] 2 T , s.t., G [f ] p (g [f ] p ) = diag{g [f ] p } = [g [f ] p ] 1 0 0 [g [f ] p ] 2 . ( 4.25)
We consider the following derivative

∂l [f ] (g [f ] ) ∂[g [f ] p ] 1 = M q=1 q>p 1 τ [f ] pq λ [f ] pq + M q=1 q<p 1 τ [f ] qp λ[f] qp (4.26)
where g

[f ] = [g [f ] T 1 , . . . , g [f ] T M ] T . Calculus and notations are specied in Appendix D, leading to [ĝ [f ] p ] 1 = a [f ] p b [f ] p (4.27) in which a [f ] p = -t [f ] * p -[g [f ] p ] 2 2S k=1 [w [f ] p ] 2k [ω [f ] * p ] 2k + 2V k=1 [ς [f ] p ] 2k [ [f ] * p ] 2k -[x [f ] p ] 1 + ρz T p [ B[f] T ] :,1 (4.28) and b [f ] p = 2S-1 k=0 [w [f ] p ] 2k+1 [ω [f ] * p ] 2k+1 + 2V -1 k=0 [ς [f ] p ] 2k+1 [ [f ] * p ] 2k+1 + ρ, (4.29)
in which we used the Matlab notation [•] :,k to refer to the k-th column.

Similarly, estimation of [g

[f ] p ] 2 is given by [g [f ] p ] 2 = ã[f] p b[f] p (4.30)
where

ã[f] p = -t [f ] * p -[g [f ] p ] 1 2S-1 k=0 [w [f ] p ] 2k+1 [ω [f ] * p ] 2k+1 + 2V -1 k=0 [ς [f ] p ] 2k+1 [ [f ] * p ] 2k+1 -[x [f ] p ] 2 +ρz T p [ B[f] T ] :,2 (4.31)
and

b[f] p = 2S k=1 [w [f ] p ] 2k [ω [f ] * p ] 2k + 2V k=1 [ς [f ] p ] 2k [ [f ] * p ] 2k + ρ. (4.32)
The proposed algorithm is based on the CG model in (3.2) and the ADMM procedure for a multi-frequency scenario. We refer to it as the Multi-frequency Robust Calibration Algorithm (MRCA) [START_REF] Ollier | Robust distributed calibration of radio interferometers with direction dependent distortions, accepted for publication[END_REF] and the global scheme is exposed in Algorithm 5.

Algorithm 5: MRCA Multi-frequency Robust Calibration Algorithm initialize: Remark: Let us note that after step 4 of MRCA in Algorithm 5, it is possible to rene the estimation: from ẑ, we deduce an estimation of the N × 1 vector z i,p , associated to the ionospheric phase delay along the path i-p. With (4.8) and (4.11), we derive N estimates for ϕ

{ θ[f] ← θ [f ] init } f ∈F , ẑ ← z init , {x [f ] ← x [f ] init } f ∈F , { Ω[f] ← Ω [f ] init } f ∈F ,{ τ [f ] ← τ [f ] init } f
[f 0 ]
i,p and perform an average. Finally, a new estimate for z can be deduced and this additional step can be directly incorporated into the estimation procedure.

In Figure 4.1, we represent the operation ow and signaling exchange between the fusion center and each local agent in our estimation procedure.

y [f ] , θ[f] , Ω[f] , τ [f ] f ∈F k-th local processor for frequency f k fusion center θ[f] f =f k θ[f]
l [f ] (θ [f ] ) = D i=1 M p=1 || Ĵ[f] i,p -G [f ] p (g [f ] p )H [f ] i,p Z [f ] i (ϕ [f ] i )F [f ] i (ϑ [f ] i )|| 2 F . ( 4
[f ] (θ [f ] , z, x [f ]
) [START_REF] Nocedal | Numerical optimization[END_REF], as no closed-form expression is available. More details can be found in our paper [START_REF]Robust calibration of radio interferometers in multi-frequency scenario[END_REF].

Simulation results

To highlight the robustness of the proposed MRCA and the benet of multi-frequency consensus optimization, we compare it with the traditional Gaussian noise assumption which is addressed with the ALS method. We also compare our multi-frequency distributed scheme with the so-called mono-frequency case where (per-channel) calibration is performed for each frequency separately, without considering a specic model of variation. Each electronic gain entry

[g

[f ]

p ] k is generated as a complex circular Gaussian random variable with mean one and variance tuning of such parameter is addressed in [START_REF] Bertsekas | Parallel and distributed computation: numerical methods[END_REF][START_REF] Yatawatta | Fine tuning consensus optimization for distributed radio interferometric calibration[END_REF]. Finally, we consider M = 8 receiving elements, D = 2 bright calibration sources and D = 4 weak background sources. We plot in Figure 4.2 the MSE of the real part of one given gain vector g [f 1 ] as a function of the SNR, for dierent number of frequencies F which are selected in the range 75 -125 MHz, the behavior being the same for other parameters of interest in {θ [f ] } f ∈F . To attain stability in convergence, less than 5 iterations are sucient in each loop. The behavior is the same for the so-called primal and dual residuals, whose convergence mostly depends on parameters ρ, N , Ñ

and initializations [START_REF] Yatawatta | Distributed radio interferometric calibration[END_REF][START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF].

The estimation results show better statistical performance when calibration is performed with the proposed MRCA in comparison with state-of-the-art methods, i.e., the conventional non-robust algorithm and/or the mono-frequency case.

As mentioned previously, multi-frequency robust calibration was also studied in the case of interferometric station (regime 3) with an ADMM-based distributed algorithm [START_REF]Robust calibration of radio interferometers in multi-frequency scenario[END_REF]. As in section 3.5.3, directional ionospheric phase delays are generated as linear functions of directional shifts α

[f ] i , which we aim to estimate. We still compare the MRCA R3 with the non-robust Gaussian case which amounts to solve a non-linear LS problem [14] and the mono-frequency case, described in Algorithm 4. The MSE for η

[f 1 ] 1
is shown as a function of the SNR in Figure 4.3 and illustrates that the multi-frequency robust calibration algorithm still provides the best statistical performance.

Conclusion

In radio astronomy, especially in regime 4, telescope design features and data processing impose major challenges in terms of calibration and the number of unknowns to estimate can be tremendously large. Thus, proposing computationally ecient algorithms is of the utmost importance and adopting distributed strategies is essential to reduce the global operational cost since the data is collected and stored among a set of computational agents. Due to separability over frequency, calibration is carried out independently for each agent and information is brought together by a fusion center which enforces consensus thanks to available constraints. The proposed algorithm is robust and exploits frequency variation of parameters of interest in regime 4 where propagation conditions dier for each element in the array and towards each source. Such scenario is adapted for calibration of large interferometric arrays but adaptation of the multi-frequency algorithm to regime 3 is also possible. In order to go even further in the simulations and obtain reconstructed images, we will use the software MeqTrees in the next chapter to simulate realistic radio telescopes and visibility measurements [START_REF] Noordam | The MeqTrees software system and its use for thirdgeneration calibration of radio interferometers[END_REF].

Chapter 5

On realistic data simulation

Presentation of MeqTrees

In this section, we give an overview of the software MeqTrees [START_REF] Noordam | The MeqTrees software system and its use for thirdgeneration calibration of radio interferometers[END_REF], which is constituted of packages well-suited for simulation and calibration of radio astronomical data for the next generation of radio telescopes. It implements measurement equations as in (2.10) with specic Jones terms and it simulates particular radio instruments. Thus, this software system is able to build numerical models by generating realistic visibility measurements but also to solve for the unknown parameters by providing benchmark calibration algorithms. Finally, it is also possible to create images after calibration thanks to the imaging package lwimager. To illustrate how the software proceeds, we show a screen capture in Figure 5.1. First of all, we need to consider a Measurement Set (MS) in which all the observation setup is dened (the antenna positions, the (u, v, w) points as described in section 2.1.1, the time slots, the frequency and the pointing center). As shown in Figure 5.1, we will be considering the KAT-7 instrument, spanning 12 hours with 60 seconds of integration time per data sample, leading to 15 120 visibility points in total.

We also need to dene the sky model with positions and intensities of all sources. As a rst step, we will be studying a set of sources taken from the Sydney University Molonglo Sky Survey (SUMSS), at 1445 MHz, using a spectral index of 0.7. A distinction is made between image-plane and uv-plane components. The rst refer to the direction independent eects which occur at the antennas while the latter stand for direction dependent eects, which the third-generation calibration aims to correct as well. We notice that many types of instrumental or atmospheric disturbances can be included, including the contribution of noise.

After compilation, we dene the size of data chunks being processed, which impacts the simulation speed. And we choose in which ouput column of the MS the simulated data are lled in (the MS can be viewed as a large array with many columns). This is shown in Figure 5.2.

Once simulations have runned and visibilities have been generated, several parameters can be set in the imaging options: column of visibilities to consider, name of output image le, image size, ... It can be seen from (2.59) that the dirty image is a convolution between the true sky image and the dirty beam. To remove sidelobes due to imperfect uv-coverage, a deconvolution process is thus required. We notice in Figure 5.2 that the CLEAN algorithm can be specied, candidates being Clark, Hogbom and CS Clean algorithms, along with parameters describing them.

The software MeqTrees also proposes a solver in the LS sense, which minimizes the dierence between observed and predicted visibilities, i.e., by tting the model to the data. To this end, the MS and the sky model, i.e., the specic sources to calibrate, are required and calibration settings need to be xed, among which the solvable parameters to consider and the type of output visibilities, as exposed in Figure 5.3. Indeed, output visibilities can be selected among the corrected data, the (un)corrected residuals and the predict. Let us consider calibration along the i-th source and for the antenna pair (p, q). If we note X pq the available data, i.e., the observed visibilities, and consider the following corrupted model (also called the predict)

V i,pq (θ) = G p (g p )K i,p C i K H i,q G H q (g q ) (5.1)
with only the typical propagation delay K i,p Jones matrix and the direction inpendent electronic gain G p term, then calibration amounts to minimize X pq -V i,pq (θ) w.r.t. θ which includes the 5.2. HANDLING OF MEQTREES: ACADEMIC EXAMPLES 67 per-antenna electronic gains g p and g q . The corrected and uncorrected residuals are, respectively, given by Ĝ-1 p (X pq -V i,pq ( θ)) ĜH -1 q (5.2) and X i,pq -V i,pq ( θ).

(5.3)

Finally, the corrected data are written as Ĝ-1 p X pq ĜH -1 q .

(5.4)

Considering the residuals enables to remove the contribution of known reference sources and unveil the weak unknown signals buried in noise and bright foregrounds. Thus, accurate calibration leads to better substraction of bright sources from the observed data and enables imaging of faint background sources which are hardly visible otherwise. We will pay particular attention to this kind of output visibilities.

Handling of MeqTrees: academic examples

To begin with, we start with a selected set of sources from the SUMSS, assumed to be unpolarized. Thus, only the intensity I is non-zero in the source brightness introduced in (2.10).

They are exposed in Table 5.1 and we will refer to them thanks to the two rst numbers in the name. seconds integration time per data sample, the full number of visibility points is 483 840 instead of 15 120. In Figure 5.5, we expose the true ideal image along with the clean and dirty images.

In the rest of the chapter, we stick to the KAT-7 instrument for computational issues.

To highlight the eects due to the presence of unknown electronic receiver gains, we make the dirty images with and without simulating periodically varying errors for gains G. In the following, we aim to compare the proposed algorithms, presented in the previous chapters, with the MeqTrees solver.

Comparison with the proposed algorithms

We wish to test the designed robust calibration techniques on realistic data sets simulated with MeqTrees and evaluate the results directly on the reconstructed residual images. To achieve this, the following steps are required

• Simulation of realistic corrupted data under MeqTrees

• Retrieve the generated data and perform calibration with the proposed algorihms

• Insert the output calibrated visibilities, e.g., the corrected residuals, into the MS

• Make the corresponding dirty or clean images with lwimager

Let us consider calibration source number 43 and D = 16 weak sources taken from the SUMSS, described in Table 5.2. They are also illustrated in Figure 5.8, without any noise or perturbation eects.

We assume that all sources, calibration and non-calibration ones, are corrupted by Faraday rotation matrices, noise is also added and we compare the calibration solver in MeqTrees with the previously introduced MRCA R3 . In Figure 5.9, we show the corrected residual images at 895 MHz in a small area surrounding the calibration source number 43, Let us note that ionospheric perturbation eects corrupt the data here. By multi-frequency non-robust calibration, we mean that a distributed strategy is adopted to handle multiple frequencies, e.g., with the ADMM procedure, but the algorithm is based on a Gaussian noise assumption.

As regards the calibration solver in MeqTrees, we also name it mono-frequency non-robust calibration as we do not exploit information across frequency (calibration is performed for each frequency separately) and the Gaussian noise model is typically assumed.

Conclusion

The MeqTrees software system is uniquely suited for generation of realistic radio astronomical data and also proposes benchmark calibration algorithms, adapted to new radio telescopes.

Various direction dependent and independent eects can be incorporated such as the receiver complex gains, the frequency bandpass and the primary beam of antennas. Each individual eect is viewed as a separable term in the software and can be corrected independently during calibration. Dirty and clean images can also be created in order to reveal the benet of calibration and the presence of weak sources of interest.

During my stay at Lab. J.-L. Lagrange in Nice (with Prof. A. Ferrari), work was conducted to master the essential features of the MeqTrees system (data simulation, calibration and imaging) and secure the linkage between the designed robust calbiration algorithms from the Matlab environment and this software widely used by the community of researchers in radio astronomy.

This was not straightforward and has taken time to develop Python scripts which handle the Chapter 6

Conclusion and prospects

The work described in this document deals with calibration of radio interferometers in non-Gaussian environment and multi-frequency scenario. Specically, a statistical study revealed that the classical Gaussian noise assumption was not suitable to model the contribution of outliers in visibility measurements, which are due to weak unknown background sources or radio frequency interferences. To address this problem, we considered the CG distribution, composed of texture and speckle parts. Deterministic and Bayesian approaches have been proposed but numerical simulations emphasized that the relaxed version (deterministic approach) was more exible and avoided model misspecications. Specically, this latter approach is based on the ML scheme. It enables to estimate iteratively the noise parameters, and more importantly, the Jones matrices which account for all perturbation eects along the signal propagation path.

Such matrices can be viewed as non-structured or structured during the estimation process. To reduce the computational burden, the EM and the BCD algorithms were used, leading to closedform expressions. For comparative purposes, traditional calibration algorithms as the classical Gaussian noise assumption and the robust Student's t-based method were considered.

To deal with multiple frequencies, we made use of a distributed and consensus strategy by exploiting the specic variation of parameters of interest w.r.t. frequency. Some variation models are available in the literature as it is the case for ionospheric disruptions, otherwise smoothness across frequency can be assumed using polynomials. To perform consensus optimization, the ADMM procedure was employed and the global computational load was distributed across a network of agents while the fusion center enforced consensus. All along the study, we focused on two particular regimes, i.e., regimes 3 and 4 as direction dependent distortions are taken into account and such scenarios are more representative of the new generation of radio interferometers.

The designed multi-frequency robust calibration algorithm was assessed against state-of-the art mono-frequency and/or non robust cases.

Finally, the use of the sotftware MeqTrees enabled to simulate realistic radio data sets with specied perturbation eects and to provide a benchmark calibration algorithm that has been compared to our work. During the calibration process, weak unknown sources contribute to the noise, aecting and changing its classical Gaussian nature but during imaging, they become the elements of interest we wish to reveal as they are not readily accessible. Thus, a better removal of calibration sources and recovering of background sources in residual images indicates a more accurate estimation of perturbation eects during calibration and enables to evaluate the performances of the designed algorithms.

Among the possible future prospects which naturally follow this work, we can mention the following points 76 CHAPTER 6. CONCLUSION AND PROSPECTS

• The temporal variation of parameters of interest could be included into the calibration process [START_REF] Van Der Tol | Self-calibration for the LOFAR radio astronomical array[END_REF]. Proposing independent solutions between blocks of time intervals is a straightforward extension and only requires to add a sum over time in the likelihood expression.

However, solving for fast-or slow-varying eects w.r.t. time, thus considering dierent time scales remains to be done. A similar approach to the one proposed in the multi-frequency scenario can be adopted.

• Compressive sensing approaches are used increasingly in radio astronomy imaging as they promote sparsity (mostly empty sky) and lead to good reconstruction results, providing competitive results with respect to traditional methods as CLEAN [START_REF] Repetti | Non-convex optimization for selfcalibration of direction-dependent eects in radio interferometric imaging[END_REF][START_REF] Dabbech | MORESANE: MOdel REconstruction by Synthesis-ANalysis Estimators-a sparse deconvolution algorithm for radio interferometric imaging[END_REF][START_REF] Bilen | Convex optimization approaches for blind sensor calibration using sparsity[END_REF][START_REF] Brossard | Calibration of radio interferometers using a sparse doa estimation framework[END_REF]. However, the Gaussian noise assumption is usually still considered. Thus, combining both aspects (sparsity and robustness) remains to be developed in future work.

• Optimization on matrix manifolds has developed signicantly recently and yields better performing algorithms when the geometric structure is exploited [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF]. Calibration is essentially an optimization problem and reformulation of the problem on a Riemannian manifold has already been proposed for radio interferometric calibration and leads to faster convergence with reduced computational cost compared to traditional calibration on Euclidean space [START_REF] Yatawatta | Radio interferometric calibration using a riemannian manifold[END_REF]. Therefore, the designed robust algorithms in this work could benet from adaptation to Riemannian manifold, using the RTR method for instance [START_REF] Absil | Trust-region methods on Riemannian manifolds[END_REF].

• The CRB is of crucial importance in estimation theory and statistics as it provides the minimal variance of unbiased estimators [START_REF] Stoica | Spectral analysis of signals[END_REF][START_REF] Kay | Fundamentals of statistical signal processing: estimation theory[END_REF]. But it can also be considered for design conguration issues in antennas of next stations [159161]. Indeed, the patterning can aect the performances of deconvolution and calibration processes and studying performance bounds can give an idea about the optimal design of stations as regards the number of antennas to consider, their positions (randomized or following a specic geometry) and the array aperture. Furthermore, studying the CRB is of particular interest, notably in order to investigate the impact of misspecied models (noise distribution, frequency variation of parameters, antenna positions) [START_REF] Mennad | Slepian-bangs-type formulas and the related misspecied Cramér-Rao Bounds for Complex Elliptically Symmetric distributions[END_REF].

• Testing the designed algorithms on real data is worth exploring and requires adaptation to very large scales [START_REF] Jongerius | An end-to-end computing model for the Square Kilometre Array[END_REF][START_REF] Wijnholds | Signal processing challenges for radio astronomical arrays[END_REF]. Indeed, with a tremendous number of radio antennas and visible radio sources in their eld-of-view, the computational cost can become a restricting factor. Eorts have been made in this work to propose algorithms with low computational complexity but deeper cost-eciency analysis would facilitate this future test phase.

Σ i,q =     α i,q β i,q 0 0 0 0 α i,q β i,q γ i,q ρ i,q 0 0 0 0 γ i,q ρ i,q

    (B.2)
in which α i,q = q * i 1 c i 1 + q * i 2 c i 3 , β i,q = q * i 1 c i 2 + q * i 2 c i 4 , γ i,q = q * i 3 c i 1 + q * i 4 c i 3 and ρ i,q = q * i 3 c i 2 + q * i 4 c i 4 .

We also obtain (3.42) where Υ i,q =     λ i,q µ i,q 0 0 ν i,q ξ i,q 0 0 0 0 λ i,q µ i,q 0 0 ν i,q ξ i,q    

(B.3)
in which λ i,q = q i 1 c i 1 + q i 2 c i 2 , µ i,q = q i 1 c i 3 + q i 2 c i 4 , ν i,q = q i 3 c i 1 + q i 4 c i 2 and ξ i,q = q i 3 c i 3 + q i 4 c i 4 .

Finally, the cost function in (3.40) We specify here the notations which lead to (4.24). To this end, we need to consider

∂l [f ] (θ [f ] ) ∂ϕ [f ] i,p and ∂h [f ] (θ [f ] ) ∂ϕ [f ] i,p
. Thus, minimization (4.19) w.r.t. phase delay ϕ [f ] i,p leads to

β [f ] i,p exp jϕ [f ] i,p + α [f ] i,p exp -jϕ [f ] i,p = 0 (C.1) in which • β [f ] i,p = M q=1 q>p j τ [f ] pq exp -jϕ [f ] i,q δ [f ] i,pq + M q=1 q<p j τ [f ] qp exp -jϕ [f ] i,q δ [f ] H i,qp + jx [f ] * i,p -jρz H i,p b [f ] * • α [f ] i,p = M q=1 q>p -j τ [f ] pq exp jϕ [f ] i,q δ [f ] H i,pq -M q=1 q<p j τ [f ] qp exp jϕ [f ] i,q δ [f ] i,qp -jx [f ] i,p + jρb [f ] T z i,p .
We also introduce

• δ [f ] i,pq =   -v [f ] H pq + D k=1 k =i s [f ] H k,pq    Ω [f ] -1 d [f ] i,pq • d [f ] i,pq = G [f ] * q (g [f ] q )H [f ] * i,q F [f ] i,q (ϑ [f ] i,q ) ⊗ G [f ] p (g [f ] p ) ⊗ H [f ] i,p F [f ] i,p (ϑ [f ] i,p ) c [f ] i .
Finally, we directly deduce (4.24) from (C.1).

Appendix D Estimation of the complex electronic gains in multi-frequency scenario

We describe here the notations introduced in the estimation of frequency dependent electronic gains. The expressions in (4.26) are given by

λ [f ] pq = t [f ] pq + w [f ] H pq I 2 ⊗ G [f ] * p [f ] pq (D.1)
in which

• t

[f ] pq = -v [f ] H pq m [f ] pq • m [f ] pq = Ω [f ] -1 G [f ] * q ⊗ E 1 w [f ] pq • E 1 = 1 0 0 0 . Besides, w [f ] pq = D i=1 H [f ] * i,q Z [f ] * i,q (ϕ i,q )F [f ] i,q (ϑ [f ] i,q ) ⊗ H [f ] i,p Z [f ] i,p (ϕ i,p )F [f ] i,p (ϑ [f ]
i,p ) c

[f ] i and

[f ] pq = G [f ]
q ⊗ I 2 m

[f ] pq .

Similarly, we have

λ[f] qp = t[f] qp + w[f] H qp G [f ] * p ⊗ I 2 ˜ [f ] qp (D.2) in which • t[f] qp = - w[f] H qp v [f ] qp • w[f] qp = Ω [f ] -1 E 1 ⊗ G [f ] q w [f ] qp ,
• ˜

[f ]

qp = I 2 ⊗ G [f ] q w [f ] qp .
The derivative in (4.26) can be written more compactly ∂l [f ] (g [f ] ) we obtain the following one-dimensional linear equation ∂l [f ] (g [f ] )

∂[g [f ] p ] 1 = t [f ] p + w [f ] H p I S ⊗ I 2 ⊗ G [f ] * p ω [f ] p + w[f] H p I V ⊗ G [f ] * p ⊗ I 2 ω[f]
∂[g [f ] p ] 1 =t [f ] p + [g [f ] * p ] 1 2S-1 k=0 [w [f ] * p ] 2k+1 [ω [f ] p ] 2k+1 + [g [f ] * p ] 1 2V -1 k=0 [ς [f ] * p ] 2k+1 [ [f ] p ] 2k+1 + [g [f ] * p ] 2 2S k=1 [w [f ] * p ] 2k [ω [f ] p ] 2k + [g [f ] * p ] 2 2V k=1 [ς [f ] * p ] 2k [ [f ] p ] 2k .
(D.4)

Considering (D.4) and

∂h [f ] (g [f ]
)

∂[g [f ]
p ] 1 leads to (4.27). Estimation of [g

[f ]
p ] 2 is similar, except that we need to consider the following quantities

• m

[f ] pq = Ω [f ] -1 G [f ] * q ⊗ E 2 w [f ] pq • w[f] qp = Ω [f ] -1 E 2 ⊗ G [f ] q w [f ] qp
• E 2 = 0 0 0 1 .

A 49 3. 9 (

 499 Per-baseline speckle component in the CG distribution B Use of the BCD algorithm for estimation of non-structured Jones matrices C Estimation of the ionospheric phase delays in multi-frequency scenario D Estimation of the complex electronic gains in multi-frequency scenario Bibliography xiv List of Figures 1 EQM de la partie réelle des 64 paramètres inconnus pour un RSB xé, avec D = 2, M = 8 et D = 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxviii 2 Régime de calibration 3, où F P et P A. . . . . . . . . . . . . . . . . . . . xxix 3 EQM d'un paramètre inconnu dans θ en fonction du RSB. . . . . . . . . . . . . . xxxi 4 Images résiduelles obtenues avec (a) la réalité terrain (b) l'algorithme de calibra-tion robuste proposé et (c) la calibration sous MeqTrees (mono-fréquence). . . . . xxxii 5 Restitution de l'amplitude en fonction de la fréquence. . . . . . . . . . . . . . . . xxxiii 2.1 Data acquisition model with interferometry (after [1]). . . . . . . . . . . . . . . . 9 2.2 Dierent uv-coverage for dierent instruments: the KAT-7 radio telescope (left) and the MeerKAT instrument (right). . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Comparison between an undistorted simulated sky (left) and a corrupted one (right) (after [2]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 Ionospheric calibration regimes (after [3, 4]). . . . . . . . . . . . . . . . . . . . . . 13 2.5 Interferometric system (after [5]). . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.6 Science data processing in the LOFAR system (after [6]). . . . . . . . . . . . . . . 19 3.1 Histograms for (a) {n } and (b) {n } (Gaussian case). . . . . . . . . . . . . . . 34 3.2 Histograms for (a){n} and (b) {n} (realistic case). . . . . . . . . . . . . . . . 35 3.3 Histograms with tted normal distribution for (a) {n } and (b) {n}. . . . . . 35 3.4 Quantile-quantile plots for (a) {n } and (b) {n}. . . . . . . . . . . . . . . . . 35 3.5 Empirical cumulative distribution plots for dierent distributions. . . . . . . . . . 36 3.6 (a) MSE of the real part of the rst 32 unknown parameters for a given SNR, (b) MSE vs. SNR for the real part of a given unknown parameter and the corresponding CRB, for D = 2 bright signal sources and M = 8 antennas. . . . . . . . . . . 48 3.7 Numerical stability for the second loop (a), and the rst loop (b) in Algorithm 3. 48 3.8 MSE for a given ionospheric phase delay as function of the SNR for dierent kinds of CG distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . a) MSE of the real part of the 64 unknown parameters for a given SNR, (b) MSE vs. SNR for the real part of a given unknown parameter, for D = 2, M = 8 and D = 8, leading to 128 real parameters of interest to estimate and 224 measurements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.10 (a) MSE of the real part of the 16 complex gains for a given SNR, (b) MSE of ζ 1 vs. SNR, for D = 2, M = 8, and D = 4, leading to 38 real parameters of interest to estimate and 224 measurements. . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.11 Recovered residual images when Jones matrices are (a) perfectly estimated (b) estimated with the RCMLE, (c) estimated with the Student's t noise model and (d) the Gaussian noise assumption (point-like sources in sky image). . . . . . . . 52 xv 3.12 Recovered residual images when Jones matrices are (a) perfectly estimated (b) estimated with the RCMLE, (c) estimated with the Student's t noise model and (d) the Gaussian noise assumption. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Communication between the k-th local processor and the fusion center. The three arrows in the center are performed sequentially and iteratively, which corresponds to the second loop in Algorithm 5. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 MSE of the real part of a given complex gain vs. SNR. . . . . . . . . . . . . . . . , for D = 2 bright signal sources, M = 8 antennas and D = 4 weak background sources. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Compile-time options of the MeqTrees software for conguration of the simulated data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Runtime options of the MeqTrees software for running the simulation and performing imaging. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Compile-time options of the MeqTrees software for calibration solver. . . . . . . . 5.4 Dirty (left) and clean (right) images in a corruption-free environment and a single source in the sky model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Ideal (left), clean (middle) and dirty (right) images of the diuse M31 source without any perturbations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Dirty images with (right) and without (left) considering disruptive G Jones gains. 5.7 Dirty images associated to the uncorrected residuals (left) and corrected residuals (right), solving for G Jones terms. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Uncorrupted dirty image with the D = 16 weak sources we wish to recover after calibration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Corrected residual images for the ideal case (left), the proposed MRCA R3 (middle) and the calibration solver in MeqTrees (right). . . . . . . . . . . . . . . . . . . . . 5.10 Recovered ux of one weak background source after calibration for F = 9 frequencies. 5.11 Corrected residual for (a) ideal case, (b) MRCA, (c) multi-frequency non-robust calibration and (d) mono-frequency non-robust calibration, around the weakest background source, i.e., source number 59. . . . . . . . . . . . . . . . . . . . . . 5.12 Corrected residual for (a) ideal case, (b) MRCA (c) multi-frequency non-robust calibration and (d) mono-frequency non-robust calibration, around three background source, i.e., sources number 25, 31 and 34. . . . . . . . . . . . . . . . . . xvi List of Tables 3.1 Overview of the dierent MAP-based estimators. . . . . . . . . . . . . . . . . . .

5 . 1

 51 Selection of point-like sources from the SUMMS. . . . . . . . . . . . . . . . . . . 5.2 Study of D = 16 weak sources from the SUMMS. . . . . . . . . . . . . . . . . . . 5.3 Recovered peak ux (in Jy) for sources in Figure 5.12. . . . . . . . . . . . . . . . xvii xviii Notations Acronyms • ADI: Alternating Direction Implicit.

4 Ω 5 τFigure 1 :

 451 Figure 1: EQM de la partie réelle des 64 paramètres inconnus pour un RSB xé, avec D = 2, M = 8 et D = 8.

Figure 2 :

 2 Figure 2: Régime de calibration 3, où F P et P A.

Figure 3 :

 3 Figure 3: EQM d'un paramètre inconnu dans θ en fonction du RSB.

Figure 4 :

 4 Figure 4: Images résiduelles obtenues avec (a) la réalité terrain (b) l'algorithme de calibration robuste proposé et (c) la calibration sous MeqTrees (mono-fréquence).

Figure 2 . 1 :

 21 Figure 2.1: Data acquisition model with interferometry (after [1]).

Figure 2 . 1 .

 21 The goal is to recover a faithful rendering of the original image, i.e., I(l, m). To this end, antennas are distributed over a specic layout in the array. Correlations are computed along each baseline vector between two antennas. Since the number of antennas is not innite, only a limited set of baselines are achieved and correlations are not measured for any point of the space (u, v), leading to an inaccurate and distorted recovered image with missing information. Let us note S(u, v) the sampling function representing actual positions for which correlations are computed, i.e., a mask CHAPTER 2. PROBLEM SETUP AND STATE-OF-THE-ART which equals 0 when no measure has been conducted at this particular position of the plane or 1 elsewhere. If we multiply the two-dimensional Fourier Transform F I(l, m) with S(u, v), the resulting data are the measured correlations. By Inverse Fourier Transform, we obtain an image, named the dirty image I D , given by

(2. 7 )

 7 As previously mentioned, synthesis mapping takes advantage of the Earth's rotation to generate a synthetic aperture. As the Earth rotates, the measured points in the (u, v) space follow elliptical tracks. The more data points we obtain, the more regular and dense is the sampling and the more accurate is the reconstructed image[48, p.24]. In Figure 2.2, we expose the uv-coverage, i.e., the set of visibility points in the (u, v) domain, for dierent simulated instruments: the KAT-7 radio telescope composed of 7 dishes on the left hand side for a total measurement time of 12 hours, and the MeerKAT consisting of 64 dishes on the right hand side during 4 hours of observation. In both cases, the integration time per visibility point is 60 seconds and simulations were performed with the Common Astronomy Software Applications (CASA) package [4951].

Figure 2 . 2 :

 22 Figure 2.2: Dierent uv-coverage for dierent instruments: the KAT-7 radio telescope (left) and the MeerKAT instrument (right).

Figure 2 . 3 :

 23 Figure 2.3: Comparison between an undistorted simulated sky (left) and a corrupted one (right)(after[START_REF] Yatawatta | Robust radio interferometric calibration[END_REF]).

Figure 2 . 4 :

 24 Figure 2.4: Ionospheric calibration regimes (after [3, 4]).

Figure 2 . 6 :

 26 Figure 2.6: Science data processing in the LOFAR system (after [6]).

  [START_REF] Noordam | The measurement equation of a generic radio telescope[END_REF] in which λ is the wavelength, L = [l 1 , . . . , l M ] T and D = [d 1 , . . . , d D ]. If the narrowband assumption holds true, geometric time delays account for phase shifts[START_REF] Ottersten | Exact and large sample maximum likelihood techniques for parameter estimation and detection in array processing[END_REF]. At each time instant, the signals measured by all antennas are collected in r(t) = As(t) + n(t)(2.25)

34 CHAPTER 3 .

 343 ROBUST CALIBRATION METHOD IN NON-GAUSSIAN ENVIRONMENT 3.1.1 Limits of the Gaussian noise modeling

Figure 3 . 1 :

 31 Figure 3.1: Histograms for (a) {n } and (b) {n } (Gaussian case).

Figure 3 . 2 :

 32 Figure 3.2: Histograms for (a) {n} and (b) {n} (realistic case).

Figure 3 . 3 :

 33 Figure 3.3: Histograms with tted normal distribution for (a) {n } and (b) {n}.

1 Figure 3 . 4 :

 134 Figure 3.4: Quantile-quantile plots for (a) {n } and (b) {n}.

Figure 3 . 5 :

 35 Figure 3.5: Empirical cumulative distribution plots for dierent distributions.

T

  and n = n T 12 , n T 13 , . . . , n T (M -1)M T includes Gaussian background noise but also the eventual presence of outliers.

pq 1 τ

 1 pq u H pq (θ)Ω -1 u pq (θ) + ln |πτ pq Ω| (3.12) 3.3. ESTIMATION OF THE TEXTURE AND SPECKLE PARAMETERS 39 in which τ = [τ 12 , τ 13 , . . . , τ (M -1)M ] T and u pq (θ) = v pq -D i=1 s i,pq (θ).

  , and the following equation to solve for the shape -BΨ(â) + pq ln(τ pq ) -B ln(b) = 0 (3.16)

(3. 18 )Algorithm 1 : 1 Obtain θ, see section 3. 4 2

 18114 IMAPE Iterative MAP estimator / K-distributed noise input : D, M , B, {C i } i∈{1,...,D} , y output : θ initialize: Ω ← Ω init , τ ← τ init , θ ← θ init while stop criterion unreached do Obtain â from solving -BΨ(a) + pq ln(τ pq ) -B ln pq τ pq Ba requires to solve -BΨ(a)pq ln(τ pq ) + B ln(b) = 0.

with D i=1 β i = 1

 1 and we note Ψ the 4B × 4B covariance matrix of n. With (3.2) and (3.3), we can write n pq ∼ CN (0, τ pq Ω).

(3. 35 )

 35 The complete data vector w = [w T 1 , . . . , w T D ] T is dened by the following 4DB × 4DB covariance matrix Ξ = bdiag{β 1 Ψ, . . . , β D Ψ}.

(3. 36 )

 36 Considering (2.45) for each source i ∈ {1, . . . , D} and [101, p. 36], the conditional expectation reads ŵi = E w i |y; θ, τ , Ω = s i (θ i ) + β i y -D l=1 s l (θ l ) .

2 M 4 Obtain

 24 -step: θi obtained as follows, i ∈ {1, . . . , D} while stop criterion unreached (referred to as the third loop) do 3 θi,p obtained from (3.43), p ∈ {1, . . . , M } Ω from (3.29) and (3.30) , 5 Obtain τ from (3.27)

Algorithm 4 : 3 Obtain

 43 RCMLE Structured case input : D, M , B, C i , β i , y, Ĵi,p as output of Algorithm 3, i ∈ {1, . . . , D} and p ∈ {1, . . . , M } output : ε initialize: ε ← ε init while stop criterion unreached do 1 Obtain θi from (3.52), i ∈ {1, . . . , D} 2 Obtain ĝp from (3.48), p ∈ {1, . . . , M } φi from (3.51), i ∈ {1, . . . , D}

Figure 3 .

 3 Figure 3.6(a), we plot the MSE of the real part of unknown entries corresponding to the rst source, leading to 32 estimated points. The Signal-to-Noise Ratio (SNR) is xed at 15 dB and the behavior remains the same for the imaginary part or the other source. For the real part of the unknown parameter vector θ, the MSE is dened as follows

Figure 3 .

 3 Figure 3.6(b), we plot the MSE of a given unknown parameter as a function of the SNR, with the corresponding CRB deduced from (3.55). Let us note that the MSE gets close to the CRB but a slight gap remains. This could be explained by our relaxed version of the exact ML estimator.

Figure 3 . 6 :

 36 Figure 3.6: (a) MSE of the real part of the rst 32 unknown parameters for a given SNR, (b) MSE vs. SNR for the real part of a given unknown parameter and the corresponding CRB, for D = 2 bright signal sources and M = 8 antennas.

Figure 3 . 7 :

 37 Figure 3.7: Numerical stability for the second loop (a), and the rst loop (b) in Algorithm 3.

CHAPTER 3 .Figure 3 . 9 :

 339 Figure 3.9: (a) MSE of the real part of the 64 unknown parameters for a given SNR, (b) MSE vs. SNR for the real part of a given unknown parameter, for D = 2, M = 8 and D = 8, leading to 128 real parameters of interest to estimate and 224 measurements.

Figure 3 .

 3 Figure 3.10: (a) MSE of the real part of the 16 complex gains for a given SNR, (b) MSE of ζ 1 vs. SNR, for D = 2, M = 8, and D = 4, leading to 38 real parameters of interest to estimate and 224 measurements.

CHAPTER 3 .Figure 3 . 11 :

 3311 Figure 3.11: Recovered residual images when Jones matrices are (a) perfectly estimated (b) estimated with the RCMLE, (c) estimated with the Student's t noise model and (d) the Gaussian noise assumption (point-like sources in sky image).

Figure 3 . 12 :

 312 Figure 3.12: Recovered residual images when Jones matrices are (a) perfectly estimated (b) estimated with the RCMLE, (c) estimated with the Student's t noise model and (d) the Gaussian noise assumption.

56 CHAPTER 4 .

 564 ROBUST DISTRIBUTED CALIBRATION IN MULTI-FREQUENCY SCENARIO

  )] T collects frequency dependent per-receiver and per-source ionospheric eects. Furthermore, we still consider the noise model in (3.2) so we also need to estimate the B × 1 texture realizations τ [f ] = [τ

Figure 4 . 1 :

 41 Figure 4.1: Communication between the k-th local processor and the fusion center. The three arrows in the center are performed sequentially and iteratively, which corresponds to the second loop in Algorithm 5.

1 4 . 62 CHAPTER 4 . 1 }

 46241 Ionospheric phase delay and Faraday rotation are randomly generated as function of the T EC, drawn from a uniform distribution U(1 × 10 17 , 5 × 10 17 ) and expressed in m -2 . We choose N = Ñ = 6 and the Lagrangian penalty factor is xed empirically as ρ = 10. Let us note that ROBUST DISTRIBUTED CALIBRATION IN MULTI-FREQUENCY SCENARIO MRCA multi-frequency robust, F = 10 MRCA multi-frequency robust, F = 3 MRCA mono-frequency robust, F = 1 ALS multi-frequency non-robust, F = 10 ALS mono-frequency non-robust, F = 1

Figure 4 . 2 :

 42 Figure 4.2: MSE of the real part of a given complex gain vs. SNR.

Figure 4 . 3 :

 43 Figure 4.3: MSE of η [f 1 ] 1 vs. SNR, for D = 2 bright signal sources, M = 8 antennas and D = 4 weak background sources.

66 CHAPTER 5 .

 665 Figure 5.1: Compile-time options of the MeqTrees software for conguration of the simulated data.

Figure 5 . 2 :

 52 Figure 5.2: Runtime options of the MeqTrees software for running the simulation and performing imaging.

Figure 5 . 3 :

 53 Figure 5.3: Compile-time options of the MeqTrees software for calibration solver.

Figure 5 . 4 :

 54 Figure 5.4: Dirty (left) and clean (right) images in a corruption-free environment and a single source in the sky model.

  Visual results are exposed in Figure5.[START_REF] Jongerius | An end-to-end computing model for the Square Kilometre Array[END_REF], where the bright source number 43 and the weaker source number 38 are considered. The position of the latter is denoted thanks to a red square around it, in the uncorrupted case. We notice that source number 38 is hardly visible but becomes fully hidden when electronic gains are included in the generated data. However, as mentioned before, it is possible to reveal the faint background by performing accurate estimation of perturbations and computing residuals. For that purpose, we still consider the previously introduced G Jones terms and perform calibration along source number 43 with the solver in MeqTrees. Dirty images associated to uncorrected and corrected residuals are shown in Figure5.7. As expected, recovering of weak source number 38 is possible.

Figure 5 . 5 :

 55 Figure 5.5: Ideal (left), clean (middle) and dirty (right) images of the diuse M31 source without any perturbations.

Figure 5 . 6 :

 56 Figure 5.6: Dirty images with (right) and without (left) considering disruptive G Jones gains.

Figure 5 . 7 :

 57 Figure 5.7: Dirty images associated to the uncorrected residuals (left) and corrected residuals (right), solving for G Jones terms.

  whose position is indicated by a red cross. We notice better removal of the calibration source for the proposed multi-frequency robust calibration algorithm while weak sources number 47 and 49 are also visible on the left hand side of the zoomed image. To illustrate how well a weak background source is recovered, we plot in Figure 5.10 the peak ux of weak source number 41 (non-calibration background source) as a function of frequency. Better ux recovering is attained with the MRCA R3 .

Figure 5 . 8 :

 58 Figure 5.8: Uncorrupted dirty image with the D = 16 weak sources we wish to recover after calibration.

Figure 5 . 9 :Figure 5 . 10 :

 59510 Figure 5.9: Corrected residual images for the ideal case (left), the proposed MRCA R3 (middle)

and 5 .

 5 [START_REF]Multisource self-calibration for sensor arrays[END_REF], we still represent the corrected residual images in restricted areas surrounding specic background sources, once calibration and subtraction of the bright source have been performed.In Figure5.11, the peak ux intensity of the recovered source (number 59) is respectively given by (a) true 0.11856 Jy, (b) 0.115736 Jy, (c) 0.105937 Jy and (d) 0.08049 Jy. The corresponding recovered peak ux of the three sources in Figure5.12 is given in Table5.3.

Figure 5 . 11 :Figure 5 . 12 :

 511512 Figure 5.11: Corrected residual for (a) ideal case, (b) MRCA, (c) multi-frequency non-robust calibration and (d) mono-frequency non-robust calibration, around the weakest background source, i.e., source number 59.

  s i,pM (θ i,p ) i θ i,p (B.[START_REF] Hamaker | Understanding radio polarimetry[END_REF] whereΣ i = [Σ T i,p+1 , • • • , Σ T i,M ] T . Likewise, we use (3.42) in si,p (θ i,p ) =    s * i,1p (θ i,p ) which Υ i = [Υ * T i,1 , • • • , Υ * T i,p-1 ] T . Inserting (B.5) and (B.6) into (B.4) and taking the derivative w.r.t. θ i,p leads to the expressions in (3.43), using the fact that A i,p and Ãi,p are Hermitian. Appendix C Estimation of the ionospheric phase delays in multi-frequency scenario

  1)p ] T . If we are considering the p-th array element, we note S = (M -p) and V = (p -1). Let us dene the permutation matrix P = [e 1 , e 3 , e 2 , e 4 ] where the 4 × 1 vector e k has zeros except at the k-th position which is equal to unity. With ς [f ] p = I V ⊗ P T w[f]

  dont les stations seront déployées en Australie et en Afrique du Sud. En tant que réel instrument d'envergure, le SKA sera constitué de millions d'antennes formant l'équivalent

d'une surface de réception d'un kilomètre carré et capables de capter des ondes radio sur un large spectre de fréquences. En sondant très rapidement le ciel avec une précision sans précédent, ce système permettra d'obtenir une qualité d'image supérieure aux installations actuelles et complètera nos connaissances sur les lois physiques régissant l'univers. Néanmoins, un tel potentiel nécessite de relever un certain nombre de dés scientiques et technologiques [911]. L'imagerie du ciel requiert notamment une estimation précise des diérentes perturbations introduites le long du trajet du signal, qu'elles soient d'origine environnementale ou instrumentale, et ce, an de reconstruire des images haute résolution, sans déformations.

  En radioastronomie, le bruit ne peut tout simplement pas être considéré comme Gaussien du fait de la présence d'interférences et de sources d'intérêt de très faible intensité, autres que les D sources de calibration supposées connues dans le modèle du ciel. Nous traitons ici la notion de

	robustesse à l'aide d'un modèle Gaussien composé pour le bruit, déni comme suit pour chaque
	paire d'antennes

T et n = n T 12 , n T 13 , . . . , n T (M -1)M T . xxiv Méthode de calibration robuste

  Estimation de θ : Conditionnellement à τ et Ω, l'estimation de θ revient à eectuer Nous distinguons alors deux cas de gures : dans le premier, les matrices de Jones sont dites non structurées et nous cherchons alors à estimer les quatre éléments complexes de chaque matrice J i,p (θ) pour i ∈ {1, . . . , D} et p ∈ {1, . . . , M }. Dans l'autre cas dit structuré, chaque matrice de Jones admet une décomposition spécique selon les perturbations physiques considérées et θ fait référence aux éléments introduits dans le modèle paramétrique[2426].

	3) l'opération suivante
	θ = argmin θ	pq	1 τ pq	u H		(12)
						)
	où t fait référence à la t-ème itération. Finalement, la contrainte spéciée préalablement sur la
	matrice de speckle conduit à la normalisation suivante
			Ωt+1 =	tr	Ωt+1 Ωt+1 .	(11)
					xxv

pq (θ)Ω -1 u pq (θ) .

  observation area is selected thanks to phase shifts introduced electronically, i.e., weighting by complex exponential terms. These antennas are called software telescopes: a beam is steered towards the direction of interest once signals captured by each individual antenna are correlated.Due to the huge number of available antennas in the array, a large collecting area is achieved and various dierent baselines are obtained, leading to better sensitivity, resolution, and a high quality image. An example of phased array telescope systems is the LOFAR[START_REF] Van Haarlem | LOFAR: The LOw-Frequency ARray[END_REF][START_REF] De Vos | The LOFAR telescope: System architecture and signal processing[END_REF][START_REF] Zarka | NenUFAR: Instrument description and science case[END_REF] whose construction started in 2006 and was completed by the Netherlands Institute for Radio Astronomy (ASTRON).Regular observations have been conducted since 2012 and approximately 50 000 antennas are spread out in stations all over Europe, particularly in the Netherlands: more than 40 stations are built in the Netherlands, 5 in Germany, 1 in France, 1 in Sweden and 1 in Great-Britain.A station is composed of multiple individual antennas: 96 low band antennas (30-80 MHz) and 48 high band tiles (120-240 MHz) which are then combined by beamforming. For now, the LOFAR is the largest radio telescope ever built (based on interferometry) and the most sensitive in the radio frequencies. Furthermore, multiple observations can be conducted in parallel, up to 488 beams simultaneously.

	8	CHAPTER 2. PROBLEM SETUP AND STATE-OF-THE-ART
		The collected signals are combined so as to simulate
	a single large dish whose diameter would be equivalent to the largest distance between two
	antennas of the array, called the aperture. More specically, the measured signals are correlated
	for each pair of antennas, spaced apart at a specic relative distance named baseline. The rst
	radio observations based on interferometry occured in 1946 when Ryle and Vonberg used two
	dipoles. Nowadays, the most popular aperture synthesis arrays cover the Westerbork Synthesis
	Radio Telescope (WSRT), built in 1970 in the Netherlands and composed of 14 dishes of 25 m
	diameter each, or the Very Large Array (VLA) in the United States inaugurated in 1980 which
	consists of 27 radio telescopes with 25 m diameter as well. For such arrays, the angular resolution
	reads	θ =	λ A	(2.2)
	where A is the aperture. Therefore, ner resolution is attainable if the distance between antennas
	is increased. However, dishes, i.e., parabolic antennas here, lack exibility. In order to target a
	specic area in the sky, they need to be moved mechanically which may slow down the response
	time and bother observation of transient phenomena.
	The new generation of radio interferometers is based on multiple exible small antennas,
	omnidirectional, of low cost and with wide eld-of-view. Contrary to mechanically based systems,

an

  1 ), . . . , r(t T )], T the number of samples and θ refers to all the unknown parameters in the model. The ALS consists in minimizing alternatingly the given cost function w.r.t. each unknown parameter, while xing the others at their previous estimates, starting from a rened enough initialization. The choice of initializations is important in the ALS as convergence to the global optimum is not ensured. However, decrease of the cost function at each iteration step can be guaranteed. As for the WALS, it consists in using an appropriate weighting term Ω s.t. the cost function becomes: ||Ω H R -R(θ) Ω|| 2 F . Better statistical performance can be achieved with the WALS, at the expense of slight increase in the computational cost.Depending on the calibration regime under study, the data model in(2.26) is dierent as well as the unknown parameters to estimate. Therefore, in what follows, we adapt the data model in (2.26) to the calibration regimes introduced in section 2.2.3

  +n . Generally, D D and only Jones matrices related to calibration sources, i.e., θ, are estimated. Parameters θ can be deduced through interpolation for instance. Gamma(ν/2, ν/2) with ν the number of degrees of freedom. The resulting marginal distribution is a Student's t dened by the following probability density function

	By removing the contribution of D strong sources, the D weakest contributions are then revealed
	if calibration has been performed accurately. One of the goals in calibration is to preserve the
	ux of the weakest sources in the resulting residual image [111].

Let us consider a Gaussian distribution N (µ, σ 2 n /γ) where µ is the mean and the variance is given by a ratio. The denominator is a random variable following a Gamma distribution, CHAPTER 2. PROBLEM SETUP AND STATE-OF-THE-ART s.t., γ ∼

Table 3 .

 3 1: Overview of the dierent MAP-based estimators.

	Texture distribution	Expression of texture pdf	Noise pdf

Table 3 .

 3

	2: Mean dierence between
	true and estimated positions for the
	D background sources.
		x	y
	RCMLE	0.0030	0.0709
	Student's t case	0.0057	0.0717
	Gaussian case	0.0066	0.0721
	Recovered intensities at the calibration and non-calibration positions are represented thanks to
	colorbars in Figure

  t+1 performed locally by each agent(4.21) where t is the iteration counter. Minimization in (4.20) requires access to local solutions from all agents at every frequency and admits the following closed-form expression

Table 5 .

 5 1: Selection of point-like sources from the SUMMS.

	Name	RA (hours h, minutes m, seconds s)	DEC (degrees • , minutes ', seconds )	Intensity (Jansky Jy)
	38_J1845M 76 19 12 16.78	-75 34 28.10	0.1097
	43_J1930M 76 19 15 54.44	-74 39 36.90	1.142

In Figure

5

.4, we expose the benet of deconvolution by showing the dirty (left) and the clean (right) images, with the single bright source number 43 in the middle which is free of any

Table 5 .

 5 3: Recovered peak ux (in Jy) for sources in Figure5.12.

	Name	Source number 31	Source number 34	Source number 25
	(a) True	0.167113	0.108398	0.122119
	(b)	0.165512	0.107406	0.122114
	(c)	0.121019	0.078651	0.0897499
	(d)	0.116302	0.0781593	0.0815132

  can be written as φ i (θ i,p ) = w i,p -s i,p (θ i,p )H A i,p w i,p -s i,p (θ i,p ) + wi,psi,p (θ i,p )H Ãi,p wi,psi,p (θ i,p ) + Constant USE OF THE BCD ALGORITHM FOR ESTIMATION OF NON-STRUCTURED JONES MATRICES• s i,p (θ i,p ) = [s T i,p(p+1) (θ i,p ), . . . , s T i,pM (θ i,p )] T • A i,p = bdiag{β i τ p(p+1) Ω, . . . , β i τ pM Ω} -1 .Furthermore, we have• wi,p = [w * T i,1p , . . . , w * T i,(p-1)p ] T • si,p (θ i,p ) = [s * T i,1p (θ i,p ), . . . , s * T i,(p-1)p (θ i,p )] T • Ãi,p = bdiag{β i τ 1p Ω * , . . . , β i τ (p-1)p Ω * } -1 .We make use of(3.41) in what follows s i,p (θ i,p ) =

	 s i,p(p+1) (θ i,p )
	 	. . .
		(B.4)

where

• w i,p = [w T i,p(p+1) , . . . , w T i,pM ] T 80 APPENDIX B.

MS data sets, to study the robustness aspects by uncluding the CG-based approaches and extend to the multi-frequency scenario. The exibility of the software also makes it possible to include our own simulated perturbation eects on the basic uncorrupted data, generated by the interferometer under study. In brief, the software MeqTrees enables to understand better what an interferometer would measure in practice but also to consider the practical feasibility of the proposed techniques in this document.

Remerciements

and include in particular the electronic gains in (2.17), the ionospheric phase shifts in (2.16) and the Faraday rotation in (2.18), which are all randomly generated here.

Bayesian vs. deterministic

We introduced dierent kinds of CG-based noise models obtained with dierent texture priors in section 3.3.1. To compare their estimation performances with the RCMLE and see how estimation of parameters of interest θ is aected, let us consider a realistic scenario with D = 2 bright calibration sources, M = 8 antennas in the sensor array and D = 4 weak background sources, resulting in non-Gaussian environment.

In Figure 3.8, we plot the MSE of the previously introduced estimators for one representative ionospheric phase delay, i.e., ϕ 1,3 here (rst source, third antenna), as a function of the SNR. To assess the robustness of the proposed methods, we also plot the MSE for traditional calibration based on a Gaussian noise assumption and the minimization of a LS cost function. As expected, this restricted hypothesis leads to poor performance while estimation results dier depending on which p(τ pq ; ϕ) is considered. Still, we notice that the lowest MSE is achieved with the RCMLE which enables to reach robustness due to its exibility and lack of misspecications, resulting in more accuracy. Therefore, the relaxed version of the ML estimator, with texture components considered as deterministic in the estimation process, will be favoured in the rest of the study.

State-of-the-art vs. deterministic

To evaluate and compare the performances, we consider the RCMLE from Algorithm 3, the robust approach based on the specic Student's t distribution [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF] and the traditional Gaussian noise assumption [14] which both compose state-of-the-art. As a rst step, we still study the non-structured case in which θ includes real and imaginary parts of all Jones matrices along the calibration source paths. The approach in [START_REF] Yatawatta | Robust radio interferometric calibration[END_REF] specically considers a Student's t noise modeling with i.i.d. entries and uses the ECME algorithm as described in section 2.4.3.3 [START_REF] Liu | ML estimation of the t distribution using EM and its extensions[END_REF][START_REF] Li | A t-distribution based particle lter for target tracking[END_REF]. Traditional calibration is based on the zero-mean white Gaussian noise assumption and solves a LS problem as in section 2.4.2 [14]. Results are plotted in Figure 3.9 for similar computation times.

Better accuracy is obtained with the RCMLE as no specic noise distribution has been chosen

and CG-based models include a wide range of dierent distributions, leading to more exibility Chapter 4 Robust distributed calibration in multi-frequency scenario

Motivation for multi-frequency calibration

The new generation of radio interferometers requires computationally ecient algorithms wich perform calibration accurately while being robust. Robustness was addressed in Chapter 3

by proposing the CG model for the noise contribution. Unknown parameters of interest were estimated for one given frequency but in practice, radio interferometers operate on a broad spectral range, resulting in large multi-frequency data from which information is retrieved. Perturbation eects, modeled by Jones matrices, as well as the noise contribution do not remain constant across frequency and their variation can be exploited during calibration in a multi-frequency scenario [START_REF] Brossard | Parallel multi-wavelength calibration algorithm for radio astronomical arrays[END_REF]. In this chapter, we focus on the frequency dependence (in addition, one might take into account the temporal variation using the same approach by considering independent solutions between time intervals). We extend the structured case, whose model was exposed in section 3.2.2 and the corresponding estimator in section 3.4.2, to regime 4 which is more adapted to large interferometric arrays. Thus, propagation conditions are variable per station and per source due to large eld-of-view and long baselines, leading to more unknowns to estimate than in regime 3 where the same part of the ionosphere is seen by all antennas in the station. Let us note that adaptation to regime 3 is straightforwardly obtained and multi-frequency calibration for the non-structured case is discussed in [START_REF] Yatawatta | Distributed radio interferometric calibration[END_REF].

Considering multiple frequency channels can signicantly increase the generated data volumes so the ow of information needs to be handled in the most optimal way. To this end, distributed calibration can be exploited with a network of computational agents. Each of them achieves calibration locally for a given subband and a fusion center collects data from all different frequencies. This method benets from data parallelism and distributed computing load, contrary to standard arrangement where a single computational agent directly operates on the full observing bandwidth and considers one frequency bin at a time [START_REF] Wijnholds | Calibration challenges for future radio telescopes[END_REF][START_REF]Multisource self-calibration for sensor arrays[END_REF]. Furthermore, making use of distributed procedures is consistent with the new radio synthesis arrays since the global bandwidth is divided into multiple frequency channels. Exploiting frequency diversity is not new and has already been tackled with bandpass solutions [START_REF] Trott | Spectral calibration requirements of radio interferometers for Epoch of Reionisation science with the SKA[END_REF] or smooth polynomials [START_REF] Yatawatta | Distributed radio interferometric calibration[END_REF] in the non-structured scenario. But the purpose of this chapter is to combine robustness with estimation of physical parameters taking into account their structured variation across frequency in regime 4.

Appendix A

Per-baseline speckle component in the

CG distribution

We describe here the corresponding expressions of (3.27) and (3.29) when a dierent speckle matrix Ω pq is assumed for each baseline p < q, p, q ∈ {1, . . . , M } 2 . In this case, the log-likelihood function is written as ln p(y|θ, τ , Ω 12 , Ω 13 , . . . ,

For each antenna pair, the texture estimate reads

while the per-baseline speckle component is given by

The remainder of the algorithm is straightforwardly obtained using (A.3).

Appendix B

Use of the BCD algorithm for estimation of non-structured Jones matrices

We present here the main steps leading to (3.43). Firstly, for sake of clarity, let us denote

T to refer to the four entries of the vectorization of source coherency matrix C i . Likewise, for the i-th source, p-th antenna and q-th antenna, respectively, we write

Using these latter notation, we obtain Title: Contributions to robust calibration methods in radio astronomy Keywords: Estimation, robust calibration, distributed optimization, radio astronomy Abstract: Accurate calibration is of critical importance for new advanced interferometric systems in radio astronomy in order to recover high resolution images with no distortions. This process consists in correcting for all environmental and instrumental effects which corrupt the observations. Most state-of-the-art calibration approaches assume a Gaussian noise model and operate mostly in an iterative manner for a mono-frequency scenario. However, in practice, the Gaussian classical noise assumption is not valid as radio frequency interference affects the measurements and multiple unknown weak sources appear within the wide field-of-view. Furthermore, considering one frequency bin at a time with a single centralized agent processing all data leads to suboptimality and computational limitations.

The goal of this thesis is to explore robustness of calibration algorithms w.r.t. the presence of outliers in a multi-frequency scenario. To this end, we propose the use of an appropriate noise model, namely, the so-called compound-Gaussian which encompasses a broad range of different heavy-tailed distributions. To combine limited computational complexity and quality of calibration, we designed an iterative calibration algorithm based on the maximum likelihood estimator under the compound-Gaussian modeling. In addition, a computationally efficient way to handle multiple subfrequency bands is to apply distributed and decentralized strategies. Thus, the global operational load is distributed over a network of computational agents and calibration amounts to solve a global constrained problem thanks to available variation models or by assuming smoothness across frequency.