Robust Calibration of Radio Interferometers in Non-Gaussian Environment
Virginie Ollier, Mohammed Nabil El Korso, Rémy Boyer, Pascal Larzabal, Marius Pesavento
aa r X i v : . [ s t a t . A P ] J u l Robust Calibration of Radio Interferometers inNon-Gaussian Environment
Virginie Ollier, Mohammed Nabil El Korso, R´emy Boyer,
Senior Member, IEEE , Pascal Larzabal,
Member, IEEE and Marius Pesavento,
Member, IEEE
Abstract —The development of new phased array systems inradio astronomy, as the low frequency array (LOFAR) and thesquare kilometre array (SKA), formed of a large number ofsmall and flexible elementary antennas, has led to significantchallenges. Among them, model calibration is a crucial stepin order to provide accurate and thus meaningful images andrequires the estimation of all the perturbation effects introducedalong the signal propagation path, for a specific source directionand antenna position. Usually, it is common to perform model cal-ibration using the a priori knowledge regarding a small numberof known strong calibrator sources but under the assumption ofGaussianity of the noise. Nevertheless, observations in the contextof radio astronomy are known to be affected by the presenceof outliers which are due to several causes, e.g., weak non-calibrator sources or man made radio frequency interferences.Consequently, the classical Gaussian noise assumption is violatedleading to severe degradation in performances. In order totake into account the outlier effects, we assume that the noisefollows a spherically invariant random distribution. Based onthis modeling, a robust calibration algorithm is presented in thispaper. More precisely, this new scheme is based on the design ofan iterative relaxed concentrated maximum likelihood estimationprocedure which allows to obtain closed-form expressions for theunknown parameters with a reasonable computational cost. Thisis of importance as the number of estimated parameters dependson the number of antenna elements, which is large for thenew generation of radio interferometers. Numerical simulationsreveal that the proposed algorithm outperforms the state-of-the-art calibration techniques.
Index Terms —Calibration, robustness, spherically invariantrandom process, relaxed concentrated maximum likelihood,Jones matrices, radio astronomy
I. I
NTRODUCTION
Radio astronomy aims to study radio emissions from thesky, in order to detect, identify new objects and observe knownstructures at higher resolution, in a specific electromagneticspectrum [1]. This fundamental thematic shines a new light onour universe, revealing more about its nature and history. Inorder to carry out particularly sensitive observations in a large
Virginie Ollier and Pascal Larzabal are with SATIE, UMR 8029,´Ecole Normale Sup´erieure de Cachan, Universit´e Paris-Saclay, Cachan,France (e-mail: [email protected], [email protected]). Virginie Ollier and R´emy Boyer are with L2S, UMR 8506,Universit´e Paris-Sud, Universit´e Paris-Saclay, Gif-sur-Yvette, France (e-mail: [email protected]). Mohammed Nabil El Korso is withLEME, EA 4416, Universit´e Paris-Ouest, Ville d’Avray, France (e-mail:[email protected]). Marius Pesavento is with Communication Sys-tems Group, Technische Universit¨at Darmstadt, Darmstadt, Germany (e-mail:[email protected]).This work was supported by MAGELLAN (ANR-14-CE23-0004-01), ONFIRE project (Jeunes Chercheurs GDR-ISIS) and ANR ASTRID projectMARGARITA. range of the spectrum, and to handle significant cosmologicalissues, largely distributed sensor arrays are currently beingbuilt or planned, such as the low frequency array (LOFAR) [2]and the square kilometre array (SKA) [3]. They will notablybe composed of a large number of relatively low-cost smallantennas with wide field of view, resulting in a large collectingarea and high resolution imaging. Nevertheless, to meet thetheoretical optimal performances of such next generation radiointerferometers, a plethora of signal processing challengesmust be overcome, among them, calibration, data reductionand image synthesis [4]–[7]. These aspects are intertwinedand must be dealt with to take advantage of the new advancedradio interferometers. As an example, lack of calibration hasdramatic effects in the image reconstruction by causing severedistortions. In this paper, we focus on calibration, whichinvolves the estimation of all unknown perturbation effectsand represents a cornerstone of the imaging step [8]–[10].Array calibration aspects have been tackled for a fewdecades in the array processing community leading to a varietyof calibration algorithms [11]–[13]. Such algorithms can beclassified into two different approaches depending on thepresence [14]–[16], or the absence [17]–[22], of one or morecooperative sources, named calibrator sources. In the radioastronomy context, calibration is commonly treated using thefirst approach as we have access to prior knowledge thanksto tables describing accuratly the position and flux of thebrightest sources [23].Following this methodology, the majority of proposed cal-ibration schemes in radio interferometry are least squares-based approaches. The state-of-the-art consists in the so-calledalternating least squares approach [24]–[27], which leads tostatistically efficient algorithm under a Gaussian model, sincethe least squares estimator is equivalent to the maximumlikelihood (ML) estimator in this case. On the other hand,expectation maximization (EM) [28]–[30] and EM-based al-gorithms, such as the space alternating generalized expectationmaximization algorithm [31], have been proposed in orderto enhance the convergence rate of the least squares-basedcalibration algorithms [32]. Nevertheless, the major drawbackof these schemes is the Gaussianity assumption which isnot realistic in the radio astronomy context. Specifically, thepresence of outliers has multiple causes, among which i) theradio frequency interferers, which corrupt the observations andare not always perfectly filtered in practice [33,34], ii) thepresence of unknown weak sources in the background [35],iii) the presence of some punctual events such as interferencedue to the Sun or due to strong sources in the sidelobes whichan also randomly create outliers [36]. To the best of ourknowledge, the proposed scheme in [35], represents the onlyalternative to the existing calibration algorithms based on aGaussian noise model.In [35], theoretical and experimental analyses have beenconducted in order to demonstrate that the effect of outliers inthe radio astronomy context can indeed be modeled by a non-Gaussian heavy-tailed distributed noise process. Nevertheless,the algorithm presented in [35] has its own limits, since thenoise is specifically modeled as a Student’s t with independentidentically distributed entries. To improve the robustness of thecalibration, we propose, in this paper, a new scheme basedon a broader class of distributions gathered under the so-called spherically invariant random noise modeling [37,38],which includes the Student’s t distribution. A sphericallyinvariant random vector (SIRV) is described as the productof a positive random variable, named texture, and the so-called speckle component which is Gaussian, resulting in atwo-scale compound Gaussian distribution [39]. The flexibilityof the SIRV modeling allows to consider non-Gaussian heavy-tailed distributed noise in the presence of outliers, but also toadaptively consider Gaussian noise in the extreme case whenthere are no outliers. Under the SIRV model, we estimatethe unknown parameters iteratively based on a relaxed MLestimator, leading to closed-form expressions for the noiseparameters while a block coordinate descent (BCD) algorithm[40,41] is designed to obtain the estimates of parameters ofinterest efficiently and at a low cost.Finally, it is worth mentioning that the parametric modelused in this paper to describe the perturbation effects isbased on the so-called Jones matrices [42,43]. Such formalismdescribes in a flexible way the conversion of the incidentelectric field into voltages. Indeed, along its propagation path,the signal is affected by various effects and transformationswhich correspond to matrix multiplications in the mathemati-cal Jones framework. Multiple distortion effects caused by theenvironment and/or the instruments can be easily incorporatedinto the model using an adequate parametrization of the Jonesmatrices. Such effects can represent, for example, the iono-spheric phase delay resulting in angular shifts, the atmosphericdistortions, the typical phase delay due to geometric pathlengthdifference, the voltage primary beam, the cross-leakage oralso the electronic gains [44,45]. For the above reasons anddue to its flexibility [1,32,42]–[44], we adopt this parametricmodel. We make a distinction between the non-structured andthe structured cases: in the first one, one total Jones matrixstands for all the effects along the full signal path while inthe second case, we regard each physical effect separatelythanks to individual Jones terms in a cumulative product. Thus,different corruptions are described by different kinds of Jonesmatrices. We emphasize that the proposed algorithm, entitledrelaxed concentrated ML estimator, is a generic algorithm asit is based on a non-structured Jones matrices formulation asa first step. However, it can be adapted to various regimesdescribing distinct calibration scenarios in which an array canoperate [46]. In this paper, we consider the specific exampleof the direction dependent distortion regime with a compactset of antennas, which we refer to as the 3DC regime. The array is therefore considered as a closely packed group ofantennas but the array elements have a wide field of view.This is particularly well-adapted for calibration of compactarrays, typically a LOFAR station.The rest of the paper is organized as follows: in SectionII, we present the data model in the context of radio astron-omy, first with non-structured Jones matrices and thereafter,we study an example of structured Jones matrices for the3DC calibration regime. In Section III, we give an overviewof the proposed robust ML estimator, based on sphericallyinvariant random process (SIRP) noise modeling. An efficientestimation procedure of the distortions introduced on eachsignal propagation path is derived in Section IV. Then, thealgorithm is adapted to the case of structured Jones matrices inSection V for the 3DC calibration regime. Finally, we providenumerical simulations in Section VI to assess the robustnessof the approach and draw our conclusions in Section VII.In this paper, we use the following notation: symbols ( · ) T , ( · ) ∗ , ( · ) H denote, respectively, the transpose, the complexconjugate and the Hermitian transpose. The Kronecker productis represented by ⊗ , E {·} denotes the expectation operator, bdiag {·} is the block-diagonal operator, whereas diag {·} con-verts a vector into a diagonal matrix. The trace and determinantoperators are, respectively, referred by tr {·} and | · | . Thesymbol I B represents the B × B identity matrix, vec( · ) stacksthe columns of a matrix on top of one another, || · || F is theFrobenius norm, while ||·|| denotes the l norm. Finally, ℜ {·} represents the real part and we note j the complex numberwhose square equals − .II. D ATA MODEL
A. Case of non-structured Jones matrices
Let us consider M antennas with known locations thatreceive D signals emitted by calibrator sources. Each antennais dual polarized and composed of two receptors, in orderto provide sensitivity to the two orthogonal polarization di-rections ( x, y ) of the incident electromagnetic plane wave.Consequently, the relation between the i -th source emissionand the measured voltage at the p -th antenna is given by[42,43,47] v i,p ( θ ) = J i,p ( θ ) s i (1)where s i = [ s i x , s i y ] T is the incoming signal, v i,p ( θ ) =[ v i,p x ( θ ) , v i,p y ( θ )] T is the generated voltage with one outputfor each polarization direction and J i,p ( θ ) denotes the so-called × Jones matrix, parametrized by the unknown vectorof interest θ . The Jones matrix models the array responseand all the perturbations introduced along the path from the i -th source to the p -th sensor. Since each propagation pathis particular, we can associate a different Jones matrix witheach source-antenna pair ( i, p ) , leading to a total number of DM Jones matrices. In this section, we consider the non-structured case where no specific perturbation model is usedto describe the physical mechanism behind each perturbationeffect and the unknown elements correspond to the entries ofall Jones matrices [32,48] (a structured example is given for3DC calibration regime, in Section II-B).2or each antenna pair, we compute the correlation of theoutput signals, resulting in the typical observations recordedby a radio interferometer. The correlation between voltages isgiven, in the case of noise free measurements, for the ( p, q ) antenna pair, by V pq ( θ ) = E ( D X i =1 v i,p ( θ ) ! D X i =1 v Hi,q ( θ ) !) = D X i =1 J i,p ( θ ) C i J Hi,q ( θ ) for p < q, p, q ∈ { , . . . , M } , (2)where the signals emitted by the sources are assumed uncor-related and the × matrix C i = E { s i s Hi } is known fromprior knowledge. Let us remark that autocorrelations are notconsidered as shown by the condition p < q in (2) (this is atypical situation in the radio astronomy context where radiointerferometric systems automatically flag the autocorrelations[2]).Using the property vec( ABC ) = ( C T ⊗ A )vec( B ) , werewrite (2) as a × vector ˜ v pq ( θ ) = vec (cid:16) V pq ( θ ) (cid:17) = D X i =1 u i,pq ( θ ) (3)in which u i,pq ( θ ) = (cid:0) J ∗ i,q ( θ ) ⊗ J i,p ( θ ) (cid:1) c i with c i =vec( C i ) . We stack all the noisy measurements within a fullvector x = h v T , v T , . . . , v T ( M − M i T ∈ C B × , where B = M ( M − denotes the total number of antenna pairs and v pq = ˜ v pq ( θ ) + n pq with n pq the noise sample at a specificantenna pair. Specifically, x reads x = D X i =1 u i ( θ ) + n (4)in which u i ( θ ) = h u Ti, ( θ ) , u Ti, ( θ ) , . . . , u Ti, ( M − M ( θ ) i T and n = h n T , n T , . . . , n T ( M − M i T is the full noise vectorwhich accounts for Gaussian noise, but also the presenceof outliers in our data. Therefore, the noise can no longerbe considered Gaussian and a robust calibration method isrequired. To investigate non-Gaussian noise modeling andencompass a broad range of noise distributions, we proposeto adopt the SIRP noise model [37,38]. Specifically, the noiseat each antenna pair is assumed to be generated as n pq = √ τ pq g pq , (5)where the positive real random variable τ pq is referred to astexture, whereas the complex speckle component g pq followsa zero-mean Gaussian distribution , i.e., g pq ∼ CN ( , Ω ) . (6)In order to remove scaling ambiguities, we impose tr { Ω } = 1 .Note that the choice of this constraint is arbitrary and does notaffect the estimates of interest as argued in [49]. Let us note that it is possible to consider a different covariance matrix Ω pq for each speckle component g pq in (6). In this case, the proposed algorithmrequires a few modifications and the corresponding expressions are presentedin Appendix A. A : array aperture V : station field-of-viewionosphere S : ionosphericirregularity scale Fig. 1. 3DC calibration regime, for which V ≫ S and S ≫ A . All receivingelements in the station see the same ionosphere part but, due to their widefield of view, a multitude of sources are visible and perturbations are highlydirection dependent. In this section, we adopted the non-structured Jones matricesformulation which is relevant in the radio astronomical context[32,48]. In this case, there is no need to specify the fullpropagation path, avoiding misspecification in the model.Besides, it is highly flexible and can be adapted to differentscenarios [46]. In the following, we present the directiondependent distortion regime with a compact set of antennas,named 3DC regime.
B. Specific case of the 3DC calibration regime
For a specific propagation path, from the i -th source to the p -th antenna, the global Jones matrix J i,p accounts for multipleeffects which can be described explicitly. Indeed, each globalmatrix can be decomposed into individual Jones terms whichstand for specific physical effects [43,44]. This way, insteadof estimating entries of all Jones matrices as done in thenon-structured case, we will estimate physically meaningfulparameters, thus reducing the total number of parameters toestimate. Introducing structured Jones matrices can be done inthe context of calibration scenarios [10,46]. In what follows,we target one specific commonly used calibration regime thatwe call 3DC calibration regime, described in Fig. 1, which iswell adapted for calibration of LOFAR on station level [26],and also for stations of the future SKA radio interferometer[10]. In this scenario, direction dependent distortions play asignificant role since individual elements in the array have3 wide field of view. Indeed, this implies different propaga-tion conditions towards distinct sources in the field of view.However, the array being relatively compact, made of similarelements, some effects might be the same for all antennas.In the following, we introduce a particular sequence of Jonesmatrices with specific parametrizations, in the context of 3DCcalibration regime [50] J i,p ( θ i,p ) = G p ( g p ) H i,p Z i,p ( α i ) F i ( ϑ i ) (7)for i ∈ { , . . . , D } , p ∈ { , . . . , M } and θ i,p =[ ϑ i , g Tp , α Ti ] T . We note H i,p the only assumed known matrixthanks to electromagnetic simulations in terms of antennaresponse and a priori knowledge given by calibrator sourceand antenna positions [43,44,50,51], whereas the remainingmatrices are explained in the following items. • Ionospheric effects :Propagation through the ionosphere, the outer layer of theearth’s atmosphere, introduces propagation delay on the signalwhich is affected by spatially variable fluctuations. If theseperturbation effects are not corrected for, the sources mayappear shifted from their intrinsic positions [10,52]. In thecase of a compact array, the ionospheric delay matrix is infact a scalar direction-dependent phase given by Z i,p ( α i ) = exp n jϕ i,p o I (8)in which ϕ i,p = η i u p + ζ i v p where α i = [ η i , ζ i ] T is thevector of unknown offsets resulting in a shift of the i -th sourcedirection and r p = [ u p , v p ] T is the vector of known antennaposition in units of wavelength.On top of that, passing through the ionosphere is associatedwith a rotation of the polarisation plane of each signal sourcearound the line of sight. We call it the ionospheric Faradayrotation matrix F i ( ϑ i ) and write it as F i ( ϑ i ) = (cid:20) cos( ϑ i ) − sin( ϑ i )sin( ϑ i ) cos( ϑ i ) (cid:21) (9)where ϑ i is the unknown Faraday rotation angle, assumedidentical for all antennas, since the array has a limited spatialextent [44]. • Instrumental effects :Individual antennas are described by electronic complex gainswhich appear in G p ( g p ) = diag { g p } with g p the unknownelectronic gain vector.Therefore, in this specific structured case, the physi-cal model parameters in (7) are collected in the vec-tor ε = P [ θ T , , θ T , , . . . , θ T D,M ] T where P isan appropriate rearrangement matrix such that ε =[ ϑ , . . . , ϑ D , g T , . . . , g TM , α T , . . . , α TD ] T .III. R OBUST CALIBRATION ESTIMATOR
This section is devoted to the design of a robust calibra-tion estimator based on the model (4). As it can be seenfrom (5), one has to specify the probability density function(pdf) of each texture parameter τ pq in order to obtain theexact ML estimates. Nevertheless, in pratical scenarios, suchprior knowledge is not available. Consequently, our idea isto make use of a relaxed version of the exact model, i.e., we assume deterministic but unknown texture realizations inthe estimation process [53,54]. This ensures more flexibilityin our algorithm as the texture distribution is not preciselydescribed and avoids any possible model misspecification,which is consistent with our motivation to design a broadrobust estimator w.r.t. the presence of outliers. On the otherhand, we adopt here an iterative procedure in order to reducethe computational cost. In doing so, the proposed algorithmsequentially updates each block of unknown parameters whilefixing the remaining parameters. This leads to a relaxedconcentrated ML based calibration estimator for which theexpression of the likelihood function, when independency isassumed between measurements, is written as f ( x | θ , τ , Ω ) = Y pq | πτ pq Ω | exp (cid:26) − τ pq a Hpq ( θ ) Ω − a pq ( θ ) (cid:27) , (10)where the vector composed of all texture realizations is τ = [ τ , τ , . . . , τ ( M − M ] T and a pq ( θ ) = v pq − ˜ v pq ( θ ) .Consequently, the log-likelihood function reads log f ( x | θ , τ , Ω ) = − B log π − X pq log τ pq − B log | Ω | − X pq τ pq a Hpq ( θ ) Ω − a pq ( θ ) . (11)In the following, we present the sequential updates of eachblock of unknown parameters, namely, τ , Ω and θ , followingthe methodology as in [55,56].
1) Derivation of ˆ τ pq : Taking the derivative of the log-likelihood function in (11) w.r.t. τ pq leads to the followingtexture estimate ˆ τ pq = 14 a Hpq ( θ ) Ω − a pq ( θ ) . (12)
2) Derivation of ˆ Ω : The derivative of the log-likelihoodfunction w.r.t. the element [ Ω ] k,l of the speckle covariancematrix, using classical differential properties [57, p. 2741],leads to − B tr n ˆ Ω − e k e Tl o + X pq τ pq a Hpq ( θ ) ˆ Ω − e k e Tl ˆ Ω − a pq ( θ ) = 0 (13)where the vector e k contains zeros except at the k -th positionwhich is equal to unity. The permutation property of the traceoperator enables to rewrite (13) as − B e Tl ˆ Ω − e k + X pq τ pq e Tl ˆ Ω − a pq ( θ ) a Hpq ( θ ) ˆ Ω − e k = 0 . (14)This finally leads to the following estimate of the specklecovariance matrix ˆ Ω = 1 B X pq τ pq a pq ( θ ) a Hpq ( θ ) . (15)Inserting (12) into (15), we obtain ˆ Ω h +1 = 4 B X pq a pq ( θ ) a Hpq ( θ ) a Hpq ( θ ) (cid:16) ˆ Ω h (cid:17) − a pq ( θ ) (16)4here h denotes the h -th iteration. Due to the introducedconstraint, we normalize the estimate of Ω by its trace, asfollows ˆ Ω h +1 = ˆ Ω h +1 tr n ˆ Ω h +1 o . (17)
3) Estimation of ˆ θ : For given Ω and τ , estimating ˆ θ isequivalent to the following minimization problem ˆ θ = argmin θ (X pq τ pq a Hpq ( θ ) Ω − a pq ( θ ) ) . (18)In the following, we aim to reduce the computational cost ofthe minimization procedure in (18) by use of the EM algo-rithm. For generality, we first adopt the non-structured Jonesmatrix formulation, which can also be specified depending onthe scenario, as shown in Section V.IV. E STIMATION OF ˆ θ FOR NON - STRUCTURED J ONESMATRICES
Estimating directly the entries of the Jones matrices avoidsspecifying any particular physical model, leading to a cali-bration algorithm which is less sensitive to model errors incomparison with algorithms based on the structured case. Dueto the possible large size of θ , a multi-dimensional parametersearch needs to be carried out to solve the optimizationproblem in (18) which requires significant computation time.To reduce this complexity, we make use of the EM algorithm.The essence of this algorithm relies on a proper parametervector partitioning as well as an adequate choice of the so-called complete data. As mentioned above, the parametersof interest θ represent the entries of all Jones matrices.Consequently, it is natural to consider the following partition θ = [ θ T , . . . , θ TD ] T = [ θ T , , . . . , θ T ,M , . . . , θ TD, , . . . , θ TD,M ] T , (19)for which the vector θ i,p ∈ R × is the parametrization ofthe path from the i -th calibrator source to the p -th sensor, i.e., J i,p ( θ ) = J i,p ( θ i,p ) . A. Use of the EM algorithm to solve (18)
The EM algorithm [28]–[30] enables to compute the MLestimates and reduce the computational cost, via the iterationof two steps. The first one is the E-step which reduces, in ourscenario, to the computation of the conditional expectationof the complete data given the observed data and the currentestimate of parameters [32,58]. Afterwards, the log-likelihoodfunction of this conditional distribution is maximized in theM-step. Therefore, this last step consists in an optimizationprocess which can be performed numerically, for instancewith the Levenberg-Marquardt (LM) algorithm [59]–[61], oranalytically if closed-form expressions are available. As weshow in the following, the optimization step is carried outw.r.t. to θ i ∈ C M × instead of θ ∈ C DM × . Therefore,the global multiple source estimation problem is reduced tomultiple single source sub-problems.
1) E-step : For the i -th source, we introduce the so-calledcomplete data vector w i such that x = D X i =1 w i (20)with w i = u i ( θ i ) + n i and n = P Di =1 n i , in which n i ∼ CN ( , β i Ψ ) . We have P Di =1 β i = 1 and Ψ is thecovariance matrix of n . Since n pq ∼ CN ( , τ pq Ω ) and withthe independence property, we obtain the following block-diagonal expression for ΨΨ = bdiag { τ Ω , . . . , τ ( M − M Ω } . (21)Let us note w = [ w T , . . . , w TD ] T the complete data vector,whose covariance matrix, denoted as Ξ , has the following form Ξ = bdiag { β Ψ , . . . , β D Ψ } . (22)With [62, p. 36], and after some calculus, the expression ofthe conditional expectation is given by ˆ w i = E { w i | x ; θ , τ , Ω } = u i ( θ i ) + β i x − D X l =1 u l ( θ l ) ! . (23)
2) M-step : The goal of this step is to estimate θ i . Once ˆ w i are computed for i ∈ { , . . . , D } , the estimated completedata vector ˆ w can be evaluated. The M-step is an optimizationproblem based on the following likelihood function where w i are independent f ( ˆ w | θ , τ , Ω ) =1 | π Ξ | exp (cid:26) − (cid:16) ˆ w − u ( θ ) (cid:17) H Ξ − (cid:16) ˆ w − u ( θ ) (cid:17)(cid:27) = D Y i =1 | πβ i Ψ | exp (cid:26) − (cid:16) ˆ w i − u i ( θ i ) (cid:17) H ( β i Ψ ) − (cid:16) ˆ w i − u i ( θ i ) (cid:17)(cid:27) . (24)To obtain an estimation of θ i , we need to minimize thefollowing cost function φ i ( θ i ) = (cid:16) ˆ w i − u i ( θ i ) (cid:17) H ( β i Ψ ) − (cid:16) ˆ w i − u i ( θ i ) (cid:17) . (25)To decrease even more the complexity cost of the proposedrobust calibration scheme, we use the BCD algorithm [40,41]in the M-step. Consequently, we obtain analytical solutions foreach single source sub-problems in (25), as shown below. B. Use of the BCD algorithm to minimize (25)
In (25), the optimization is performed w.r.t. θ i . However,this unknown parameter vector can be partitioned accordingto the antennas, as expressed in (19). In the following, weperform the optimization of the cost function w.r.t. each θ i,p ( p -th antenna), given the current estimates of all the other θ i,q with q = p . This leads to a closed-form expression of ˆ θ i,p as function of θ i,q for q = p and the optimization process isrepeated for each component vector θ i,p for p ∈ { , . . . , M } .If we want to minimize (25) w.r.t. block-coordinate vector θ i,p , we notice that only a subset of u i ( θ i ) is actually5ependent on θ i,p , i.e., { u i,pq } for q > p, q ∈ { , . . . , M } and { u i,qp } for q < p, q ∈ { , . . . , M } . Therefore, (25)reads φ i ( θ i,p ) = M X q =1 q>p (cid:16) w i,pq − u i,pq ( θ i,p ) (cid:17) H ( β i τ pq Ω ) − (cid:16) w i,pq − u i,pq ( θ i,p ) (cid:17) + M X q =1 q
We recall that the output of Algorithm 1 is the estimateof each Jones matrix denoted by ˆ J i,p for i ∈ { , . . . , D } and p ∈ { , . . . , M } . In the following, we consider the data modelin (7) for the specific 3DC calibration regime, and intend toestimate the unknown parameter vector of interest ε in a Algorithm 1:
Relaxed concentrated ML based calibrationalgorithm input : D , M , B , C i , β i , x output : ˆ θ initialize: ˆ Ω ← Ω init , ˆ τ ← τ init , ˆ θ ← θ init while stop criterion unreached do while stop criterion unreached do E-step: ˆ w i obtained from (23), i ∈ { , . . . , D } M-step: ˆ θ i obtained as follows, i ∈ { , . . . , D } while stop criterion unreached do ˆ θ i,p obtained from (29), p ∈ { , . . . , M } end end Obtain ˆ Ω from (16) and (17), Obtain ˆ τ from (12) end sequential manner. To do so, we use an iterative estimationprocedure by optimizing a cost function w.r.t. one parameterwhile fixing the others.
1) Estimation of g p : The diagonal elements of the gainmatrix are given by ˆ g p = argmin g p κ ( g p ) (30)where κ ( g p ) = P Di =1 || ˆ J i,p − G p ( g p ) H i,p Z i,p F i || F . Werewrite the cost function as κ ( g p ) = D X i =1 Tr n(cid:16) ˆ J i,p − G p ( g p ) R i,p (cid:17)(cid:16) ˆ J i,p − G p ( g p ) R i,p (cid:17) H o (31)in which R i,p = H i,p Z i,p F i . The derivation of κ ( g p ) w.r.t.the k -th element [ g p ] k leads to ∂κ ( g p ) ∂ [ g p ] k = D X i =1 Tr n − e k e Tk R i,p ˆ J Hi,p + e k e Tk R i,p R Hi,p G Hp o . (32)Let us denote X i,p = R i,p ˆ J Hi,p and W i,p = R i,p R Hi,p . From(32), we deduce the equation satisfied by the gain matrix D X i =1 [ X i,p ] k,k = D X i =1 [ W i,p ˆ G Hp ] k,k = D X i =1 [ W i,p ] k,k [ ˆ G ∗ p ] k,k (33)for k ∈ { , } , since G p is a diagonal matrix. Therefore, eachcomplex gain element is estimated as [ˆ g p ] k = (cid:16) D X i =1 [ W ∗ i,p ] k,k (cid:17) − D X i =1 [ X ∗ i,p ] k,k . (34)
2) Estimation of α i : We first need to estimate ϕ i,p (werecall that ϕ i,p = η i u p + ζ i v p ). This is done as follows ˆ ϕ i,p = argmin ϕ i,p ˜ κ ( ϕ i,p ) (35)6here ˜ κ ( ϕ i,p ) = || ˆ J i,p − G p H i,p Z i,p ( ϕ i,p ) F i || F . Taking thederivative of ˜ κ ( ϕ i,p ) w.r.t. ϕ i,p and setting the result to zero,we obtain Tr n j exp − j ˆ ϕ i,p ˆ J i,p F Hi H Hi,p G Hp − j exp j ˆ ϕ i,p G p H i,p F i ˆ J Hi,p o = 0 (36)which leads to exp n j ˆ ϕ i,p o = Tr n M i,p o Tr n M Hi,p o (37)where M i,p = ˆ J i,p F Hi H Hi,p G Hp .In the case of a compact array, we can write for the i -thsource ϕ Ti = ˆ α Ti Λ (38)where ϕ i = [ ˆ ϕ i, , . . . , ˆ ϕ i,M ] T and Λ = (cid:20) u , . . . , u M v , . . . , v M (cid:21) .Therefore, estimation of the directional shifts due to propaga-tion in the ionosphere is given by ˆ α Ti = ϕ Ti Λ H " P Mp =1 v p − P Mp =1 u p v p − P Mp =1 v p u p P Mp =1 u p Mp =1 u p P Mp =1 v p − ( P Mp =1 u p v p ) . (39)
3) Estimation of ϑ i : We consider the following minimiza-tion problem ˆ ϑ i = argmin ϑ i M X p =1 || ˆ J i,p − G p H i,p Z i,p F i ( ϑ i ) || F . (40)We assume a large number of antennas M while the numberof calibrator sources D is relatively reduced, such that ob-servations outnumber unknown parameters. For each source,the 1D optimization in (40) can be computed in a reasonablecomputational time through a classical data grid search or aNewton type algorithm.Finally, the proposed algorithm for the structured Jonesmatrices case regarding 3DC calibration regime is given inAlgorithm 2. Algorithm 2:
Case of structured Jones matrices input : D , M , B , C i , β i , x , ˆ J i,p as output ofAlgorithm 1, i ∈ { , . . . , D } and p ∈ { , . . . , M } output : ˆ ε initialize: ˆ ε ← ε while stop criterion unreached do Obtain ˆ ϑ i from (40), i ∈ { , . . . , D } Obtain ˆ g p from (34), p ∈ { , . . . , M } Obtain α i = [ η i , ζ i ] T from (39), i ∈ { , . . . , D } end VI. N
UMERICAL SIMULATIONS
In this part, we first aim to assess the statistical performanceof Algorithm 1, when the noise model matches our noiseassumption, i.e., a SIRP noise modeling. Afterwards, weintend to study our proposed algorithm, in a more realisticscenario where outliers are present. More specifically, the skyis composed of D known bright calibrator sources but also D ′ weak non-calibrator sources which are absorbed in thenoise component and act as outliers. Under such assumption,we compare our scheme with the recently introduced robustcalibration approach based on Student’s t [35] and with thetraditional Gaussian cases [32]. Finally, we apply Algorithm 2to the introduced 3DC calibration regime where Jones matricesare structured. A. Numerical results under SIRP noise
The unknown parameters to estimate, θ given in (19),correspond to the real and imaginary parts of the entries ofall Jones matrices. The additive noise in (4) is assumed tofollow a SIRP as in (5) and the number of Monte Carlo runsis set to 100.In order to evaluate the estimation performance of therelaxed concentrated ML based calibration algorithm, we fixthe noise distribution, e.g., as a Student’s t and make use ofthe Cram´er-Rao bound (CRB) [65]. To do so, each randomtexture component is supposed to follow an inverse gammadistribution [66]. As an example, τ pq ∼ IG ( ν/ , ν/ , (41)with ν degrees of freedom [67] and we choose for example [ Ω ] k,l = σ . | k − l | exp j π ( k − l ) .The covariance inequality principle states that, under quitegeneral/weak conditions, the variance satisfies MSE([ ˆ θ ] k ) = E n(cid:16) [ ˆ θ ] k − [ θ ] k (cid:17) o > [CRB( θ )] k,k (42)where the CRB is given as the inverse of the Fisher informa-tion matrix (FIM) F . A Slepian-Bangs type formula of the FIMfor SIRP observations is given in [68] which can be adaptedto our case and reads [ F ] k,l = 2 ν + 4 ν + 5 X pq ℜ ( ∂ ˜ v Hpq ( θ ) ∂ [ θ ] k Ω − ∂ ˜ v pq ( θ ) ∂ [ θ ] l ) . (43)Note that the noise parameters are decoupled from the param-eters of interest. Consequently, only the part corresponding tothe latter is kept in the FIM expression.First, in Fig. 2(a), we plot the mean square error (MSE)of the real part of each unknown parameter, obtained withAlgorithm 1, for a signal-to-noise ratio SNR = 15dB . Weonly plot the parameters relative to one given source, thebehavior being the same for any source. We also comparethe MSE of one given parameter, as a function of the SNR,to its corresponding CRB in Fig. 2(b). This enables to assessthe statistical performance of Algorithm 1, and we notice thatthe MSE approaches the CRB. The small gap between thebound and the algorithm is explained by the relaxed naturehypothesis used in the design of Algorithm 1. Indeed, we7
10 20 30 40 50 60 70 parameter index -7 -6 -5 -4 M SE f o r r ea l pa r t o f θ (a) -5 0 5 10 15 SNR (dB) -7 -6 -5 -4 M SE f o r ea l pa r t o f θ Proposed algorithmCRB (b)Fig. 2. (a) MSE of the real part of the first unknown parameters for agiven SNR, (b) MSE vs. SNR for the real part of a given unknown parameterand the corresponding CRB, for D = 2 bright signal sources and M = 8 antennas, leading to real unknown parameters of interest to estimate and measurements. have assumed unknown and deterministic texture parameterswhen we derived the estimates using Algorithm 1, but inthe data model, these parameters are generated as randomvariables following inverse gamma distribution and the CRBwas derived using the prior of the pdf of the texture.Second, we now aim to investigate numerically the conver-gence properties of our algorithm, which is composed of 3loops. In each of these loops, θ is updated at each iteration.We therefore consider the following quantity ǫ h ℜ{ θ } = ||ℜ (cid:8) θ h − θ h − (cid:9) || (44)where h refers to the h -th iteration. In Fig. 3 we present theconvergence rate of loops described in Algorithm 1 at line 1and 2 (the analysis of convergence of the third loop, given inline 5, has the same behavior as the loop in line 2 and thus,is not reported here). We note that around 5 iterations arerequired in loop 2 to attain convergence while approximately20 iterations are needed for the algorithm to be stable, in −5 iteration ε R e ( θ ) (a) −3 iteration ε R e ( θ ) (b)Fig. 3. ǫ h ℜ{ θ } as function of the h -th iteration, for loop in line 2 (a), and inline 1 (b) from Algorithm 1. loop 1. Nevetheless, in simulations, we notice that only 3 to4 iterations are sufficient to get close to the CRB. B. Numerical results using a realistic model
Let us investigate the robustness of our proposed calibrationprocedure in a realistic situation, and compare it with thestate-of-the-art. To do so, we consider D calibrator sources, D ′ weak outlier sources and Gaussian background noise inour data model. The real parameters of interest to estimatestill correspond to the real and imaginary entries of theJones matrices associated to the calibrator sources paths. Inthe following, we compare Algorithm 1 with i) the calibra-tion approach exposed in [35] which assumes a Student’st noise modeling using the so-called expectation-conditionalmaximization either algorithm [69,70] and ii) the traditionalcalibration scheme based on zero-mean white Gaussian noisemodeling using a least squares approach [32]. The comparativeresults are plotted in Fig. 4, for the same computation times.We notice better estimation performance with Algorithm 1,since we did not specify any particular noise distribution in8
10 20 30 40 50 60 70 parameter index -2 -1 M SE f o r r ea l pa r t o f θ SIRP assumption (proposed algorithm)Gaussian assumption (Kazemi et al., 2013 - Madsen et al., 2004)Student's t assumption (Kazemi et al., 2013) (a) -5 0 5 10 15
SNR (dB) -2 -1 M SE f o r r ea l pa r t o f θ SIRP assumption (proposed algorithm)Gaussian assumption (Kazemi et al., 2013 - Madsen et al., 2004)Student's t assumption (Kazemi et al., 2013) (b)Fig. 4. (a) MSE of the real part of the unknown parameters for a givenSNR, (b) MSE vs. SNR for the real part of a given unknown parameter, for D = 2 , M = 8 and D ′ = 8 , leading to real parameters of interest toestimate and measurements. our procedure and a SIRP includes many different types ofdistributions. Furthermore, no assumption was made aboutindependent entries in the noise vector, thus ensuring moreflexibility and robustness. C. Structured case
In order to compare our proposed global algorithm (Algo-rithm 1 followed by Algorithm 2) with the approach based onStudent’s t [35] and the Gaussian case [32] which were bothintroduced in the non-structured case, we apply Algorithm 2on the output of these two latter algorithms.Each Jones matrix is generated according to (7) with g = [ g T , . . . , g TM ] T . In all exposed simulations, we considersimilar computation times for the three presented patterns. Weshow the results only for the complex gains, cf. Fig. 5(a), andthe source offset ζ , cf. Fig 5(b), due to lack of space, thebehavior being the same for the other parameters, i.e., ϑ , ϑ , η , η and ζ . In the case of structured Jones matrices,adapted to the 3DC calibration regime, and in the presence of parameter index -4 -3 -2 -1 M SE f o r r ea l pa r t o f g SIRP assumption (proposed algorithm)Student's t assumption (Kazemi et al., 2013)Gaussian assumption (Kazemi et al., 2013 - Madsen et al., 2004) (a) -5 0 5 10 15
SNR (dB) -5 -4 -3 -2 M SE f o r ζ SIRP assumption (proposed algorithm)Student's t assumption (Kazemi et al., 2013)Gaussian assumption (Kazemi et al., 2013 - Madsen et al., 2004) (b)Fig. 5. (a) MSE of the real part of the complex gains for a given SNR,(b) MSE of ζ vs. SNR, for D = 2 , M = 8 , and D ′ = 4 , leading to realparameters of interest to estimate and measurements. outliers, we still notice the better performances of Algorithm 1,compared to the state-of-the-art. This is expected since betterestimation of the Jones entries leads to better estimation of thephysical parameters describing the structured Jones matrices.VII. C ONCLUSION
In this paper, we proposed a robust calibration techniquewhere perturbation effects are modeled thanks to Jones ma-trices. To deal with the presence of outliers in our data,the introduced ML estimation method is based on SIRPnoise modeling, leading to a relaxed concentrated ML basedcalibration algorithm. Numerical simulations show that theproposed algorithm is more robust to the presence of outliersin comparison with the state-of-the-art, for both non-structuredand structured Jones matrices, with a reasonable computationalcomplexity. A
PPENDIX
AWe describe here the corresponding expressions of (12) and(16) when we assume a different Ω pq for p < q, p, q ∈ , . . . , M } . In this case, the log-likelihood function is writtenas log f ( x | θ , τ , Ω , Ω , . . . , Ω ( M − M ) = − B log π − X pq log τ pq − X pq log | Ω pq | − X pq τ pq a Hpq ( θ ) Ω − pq a pq ( θ ) . (45)For each antenna pair, the texture estimate is given by ˆ τ pq = 14 a Hpq ( θ ) Ω − pq a pq ( θ ) (46)while the speckle covariance estimate reads ˆ Ω h +1 pq = 4 a pq ( θ ) a Hpq ( θ ) a Hpq ( θ ) (cid:16) ˆ Ω hpq (cid:17) − a pq ( θ ) . (47)The remainder of the algorithm is straightforwardly obtainedusing (47). A PPENDIX
BWe present here the steps to obtain (29). Firstly, for sake ofclarity, let us denote c i = [ c i , c i , c i , c i ] T to refer to the fourentries of the vectorization of source coherency matrix C i .Likewise, for the i -th source, we write J i,p ( θ i,p ) = (cid:20) p i p i p i p i (cid:21) for the p -th antenna and J i,q ( θ i,q ) = (cid:20) q i q i q i q i (cid:21) for the q -th antenna, i.e., θ i,p = [ p i , p i , p i , p i ] T and θ i,q =[ q i , q i , q i , q i ] T . Using these latter notation, we obtain (27)where Σ i,q = α i,q β i,q α i,q β i,q γ i,q ρ i,q γ i,q ρ i,q (48)in which α i,q = q ∗ i c i + q ∗ i c i , β i,q = q ∗ i c i + q ∗ i c i , γ i,q = q ∗ i c i + q ∗ i c i and ρ i,q = q ∗ i c i + q ∗ i c i .We also obtain (28) where Υ i,q = λ i,q µ i,q ν i,q ξ i,q λ i,q µ i,q ν i,q ξ i,q (49)in which λ i,q = q i c i + q i c i , µ i,q = q i c i + q i c i , ν i,q = q i c i + q i c i and ξ i,q = q i c i + q i c i .Finally, the cost function in (26) can be written as φ i ( θ i,p ) = (cid:16) w i,p − u i,p ( θ i,p ) (cid:17) H A i,p (cid:16) w i,p − u i,p ( θ i,p ) (cid:17) + (cid:16) ˜ w i,p − ˜ u i,p ( θ i,p ) (cid:17) H ˜ A i,p (cid:16) ˜ w i,p − ˜ u i,p ( θ i,p ) (cid:17) + Constant (50)where w i,p = [ w Ti,p ( p +1) , . . . , w Ti,pM ] T , u i,p ( θ i,p ) =[ u Ti,p ( p +1) ( θ i,p ) , . . . , u Ti,pM ( θ i,p )] T and A i,p =bdiag { β i τ p ( p +1) Ω , . . . , β i τ pM Ω } − .Furthermore, we have ˜ w i,p = [ w ∗ T i, p , . . . , w ∗ T i, ( p − p ] T , ˜ u i,p ( θ i,p ) = [ u ∗ T i, p ( θ i,p ) , . . . , u ∗ T i, ( p − p ( θ i,p )] T and ˜ A i,p =bdiag { β i τ p Ω ∗ , . . . , β i τ ( p − p Ω ∗ } − . We make use of (27) in what follows u i,p ( θ i,p ) = u i,p ( p +1) ( θ i,p ) ... u i,pM ( θ i,p ) = Σ i,p +1 θ i,p ... Σ i,M θ i,p = Σ i θ i,p (51)where Σ i = [ Σ Ti,p +1 , · · · , Σ Ti,M ] T . Likewise, we use (28) in ˜ u i,p ( θ i,p ) = u ∗ i, p ( θ i,p ) ... u ∗ i, ( p − p ( θ i,p ) = Υ ∗ i, θ i,p ... Υ ∗ i,p − θ i,p = Υ i θ i,p (52)in which Υ i = [ Υ ∗ T i, , · · · , Υ ∗ T i,p − ] T . Inserting (51) and (52)into (50) and taking the derivative w.r.t. θ i,p leads to theexpressions in (29), using the fact that A i,p and ˜ A i,p areHermitian. R EFERENCES[1] A. R. Thompson, J. M. Moran, and G. W. Swenson Jr.,
Interferometryand synthesis in radio astronomy, 2nd ed.
John Wiley & Sons, 2001.[2] M. P. Van Haarlem et al. , “LOFAR: The LOw-Frequency ARray,”
Astronomy & Astrophysics , vol. 556, no. A2, 2013.[3] P. E. Dewdney, P. J. Hall, R. T. Schilizzi, and T. J. L. Lazio, “TheSquare Kilometre Array,”
Proceedings of the IEEE , vol. 97, no. 8, pp.1482–1496, 2009.[4] S. J. Wijnholds and A.-J. van der Veen, “Fundamental imaging limitsof radio telescope arrays,”
IEEE Journal of Selected Topics in SignalProcessing , vol. 2, no. 5, pp. 613–623, 2008.[5] U. Rau, S. Bhatnagar, M. A. Voronkov, and T. J. Cornwell, “Advances incalibration and imaging techniques in radio interferometry,”
Proceedingsof the IEEE , vol. 97, no. 8, pp. 1472–1481, 2009.[6] M. de Vos, A. W. Gunst, and R. Nijboer, “The LOFAR telescope: Systemarchitecture and signal processing,”
Proceedings of the IEEE , vol. 97,no. 8, pp. 1431–1437, 2009.[7] A.-J. van der Veen and S. J. Wijnholds, “Signal processing tools for radioastronomy,” in
Handbook of Signal Processing Systems . Springer, 2013,pp. 421–463.[8] D. A. Mitchell, L. J. Greenhill, R. B. Wayth, R. J. Sault, C. J. Lonsdale,R. J. Cappallo, M. F. Morales, and S. M. Ord, “Real-time calibrationof the Murchison Widefield Array,”
IEEE Journal of Selected Topics inSignal Processing , vol. 2, no. 5, pp. 707–717, 2008.[9] S. J. Wijnholds, A.-J. van der Veen, F. De Stefani, E. La Rosa, andA. Farina, “Signal processing challenges for radio astronomical arrays,”in
International Conference on Acoustics, Speech and Signal Processing ,Florence, Italy, 2014, pp. 5382–5386.[10] S. J. Wijnholds, S. van der Tol, R. Nijboer, and A.-J. van der Veen, “Cal-ibration challenges for future radio telescopes,”
IEEE Signal ProcessingMagazine , vol. 27, no. 1, pp. 30–42, 2010.[11] D. R. Fuhrmann, “Estimation of sensor gain and phase,”
IEEE Trans-actions on Signal Processing , vol. 42, no. 1, pp. 77–87, 1994.[12] B. C. Ng and C. M. S. See, “Sensor-array calibration using a maximum-likelihood approach,”
IEEE Transactions on Antennas and Propagation ,vol. 44, no. 6, pp. 827–835, 1996.[13] B. P. Ng, J. P. Lie, M. H. Er, and A. Feng, “A practical simple geometryand gain/phase calibration technique for antenna array processing,”
IEEETransactions on Antennas and Propagation , vol. 57, no. 7, pp. 1963–1972, 2009.[14] J. T. H. Lo and S. L. Marple, “Eigenstructure methods for array sensorlocalization,” in
International Conference on Acoustics, Speech, andSignal Processing , Dallas, Texas, 1987, pp. 2260–2263.[15] B. C. Ng and W. Ser, “Array shape calibration using sources in knownlocations,” in
ICCS/ISITA.’Communications on the Move’ , Singapore,1992, pp. 836–840.[16] B. C. Ng and A. Nehorai, “Active array sensor localization,”
Signalprocessing , vol. 44, no. 3, pp. 309–327, 1995.[17] Y. Rockah and P. Schultheiss, “Array shape calibration using sourcesin unknown locations–part i: Far-field sources,”
IEEE Transactions onAcoustics, Speech, and Signal Processing , vol. 35, no. 3, pp. 286–299,1987.
18] A. J. Weiss and B. Friedlander, “Array shape calibration using sources inunknown locations-a maximum likelihood approach,”
IEEE Transactionson Acoustics, Speech, and Signal Processing , vol. 37, no. 12, pp. 1958–1966, 1989.[19] B. Friedlander and A. J. Weiss, “Direction finding in the presence ofmutual coupling,”
IEEE Transactions on Antennas and Propagation ,vol. 39, no. 3, pp. 273–284, 1991.[20] M. P. Wylie, S. Roy, and H. Messer, “Joint DOA estimation and phasecalibration of linear equispaced (LES) arrays,”
IEEE Transactions onSignal Processing , vol. 42, no. 12, pp. 3449–3459, 1994.[21] L. Qiong, G. Long, and Y. Zhongfu, “An overview of self-calibration insensor array processing,” in , Beijing, China, 2003, pp. 279–282.[22] B. P. Flanagan and K. L. Bell, “Array self-calibration with large sensorposition errors,”
Signal Processing , vol. 81, no. 10, pp. 2201–2214, 2001.[23] J. Baars, R. Genzel, I. Pauliny-Toth, and A. Witzel, “The absolutespectrum of CAS A-an accurate flux density scale and a set of secondarycalibrators,”
Astronomy and Astrophysics , vol. 61, no. 1, pp. 99–106,1977.[24] A.-J. Boonstra and A.-J. van der Veen, “Gain calibration methods forradio telescope arrays,”
IEEE Transactions on Signal Processing , vol. 51,no. 1, pp. 25–38, 2003.[25] S. J. Wijnholds and A.-J. van der Veen, “Multisource self-calibration forsensor arrays,”
IEEE Transactions on Signal Processing , vol. 57, no. 9,pp. 3512–3522, 2009.[26] S. J. Wijnholds, “Fish-eye observing with phased array radio telescopes,”Ph.D. dissertation, Technische Universiteit Delft, Delft, The Netherlands,2010.[27] S. Salvini and S. J. Wijnholds, “Fast gain calibration in radio astronomyusing alternating direction implicit methods: Analysis and applications,”
Astronomy & Astrophysics , vol. 571, p. A97, 2014.[28] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihoodfrom incomplete data via the EM algorithm,”
Journal of the RoyalStatistical Society. Series B (Methodological) , vol. 39, no. 1, pp. 1–38,1977.[29] T. K. Moon, “The Expectation-Maximization algorithm,”
IEEE SignalProcessing Magazine , vol. 13, no. 6, pp. 47–60, 1996.[30] G. J. McLachlan and T. Krishnan,
The EM algorithm and extensions,2nd ed.
John Wiley & Sons, 2008.[31] J. A. Fessler and A. O. Hero, “Space-alternating generalized expectation-maximization algorithm,”
IEEE Transactions on Signal Processing ,vol. 42, no. 10, pp. 2664–2677, 1994.[32] S. Yatawatta, S. Zaroubi, G. de Bruyn, L. Koopmans, and J. Noordam,“Radio interferometric calibration using the SAGE algorithm,” in , Marco Island, FL, 2009, pp. 150–155.[33] J. Raza, A.-J. Boonstra, and A.-J. Van der Veen, “Spatial filtering ofRF interference in radio astronomy,”
IEEE Signal Processing Letters ,vol. 9, no. 2, pp. 64–67, 2002.[34] S. Van Der Tol and A.-J. Van Der Veen, “Performance analysis of spatialfiltering of RF interference in radio astronomy,”
IEEE Transactions onSignal Processing , vol. 53, no. 3, pp. 896–910, 2005.[35] S. Yatawatta and S. Kazemi, “Robust radio interferometric calibration,”in
International Conference on Acoustics, Speech and Signal Processing ,Florence, Italy, 2014, pp. 5429–5433.[36] A.-J. Boonstra, “Radio frequency interference mitigation in radio as-tronomy,” Ph.D. dissertation, Technische Universiteit Delft, Delft, TheNetherlands, 2005.[37] E. Jay, “D´etection en environnement non gaussien,” Ph.D. dissertation,Universit´e de Cergy Pontoise, 2002.[38] K. Yao, “Spherically invariant random processes: Theory and applica-tions,” in
Communications, Information and Network Security , V. Bhar-gava et al. , Eds. Dordrecht, The Netherlands: Kluwer Academic:Springer, 2003, pp. 315–331.[39] J. Wang, A. Dogandzic, and A. Nehorai, “Maximum likelihood es-timation of compound-gaussian clutter and target parameters,”
IEEETransactions on Signal Processing , vol. 54, no. 10, pp. 3884–3898, 2006.[40] J. Friedman, T. Hastie, H. H¨ofling, R. Tibshirani et al. , “Pathwisecoordinate optimization,”
The Annals of Applied Statistics , vol. 1, no. 2,pp. 302–332, 2007.[41] M. Hong, M. Razaviyayn, Z.-Q. Luo, and J.-S. Pang, “A unifiedalgorithmic framework for block-structured optimization involving bigdata: With applications in machine learning and signal processing,”
IEEESignal Process. Mag. , vol. 33, no. 1, pp. 57–77, 2016.[42] J. P. Hamaker, J. D. Bregman, and R. J. Sault, “Understanding radiopolarimetry. I. Mathematical foundations,”
Astronomy and AstrophysicsSupplement Series , vol. 117, no. 1, pp. 137–147, 1996. [43] O. M. Smirnov, “Revisiting the radio interferometer measurement equa-tion. I. A full-sky Jones formalism,”
Astronomy & Astrophysics , vol.527, no. A106, 2011.[44] J. E. Noordam, “The measurement equation of a generic radio telescope,AIPS++ implementation note nr 185,” ASTRON, Tech. Rep., 1996.[45] J. E. Noordam and O. M. Smirnov, “The MeqTrees software systemand its use for third-generation calibration of radio interferometers,”
Astronomy & Astrophysics , vol. 524, no. A61, 2010.[46] C. Lonsdale, “Calibration approaches,” LFD memo 015, Tech. Rep.,2004.[47] L. M. Ker, “Radio AGN evolution with low frequency radio surveys,”Ph.D. dissertation, The University of Edinburgh, 2012.[48] C. D. Nunhokee, “Link between ghost artefacts, source suppressionand incomplete calibration sky models,” Ph.D. dissertation, RhodesUniversity, Rhodes, South Africa, 2015.[49] X. Zhang, M. N. El Korso, and M. Pesavento, “MIMO radar targetlocalization and performance evaluation under SIRP clutter,”
SignalProcessing Journal , vol. 130, pp. 217–232, 2017.[50] S. Yatawatta, “Reduced ambiguity calibration for LOFAR,”
Experimen-tal Astronomy , vol. 34, no. 1, pp. 89–103, 2012.[51] S. Kazemi, P. Hurley, O. ¨Oc¸al, and G. Cherubini, “Blind calibration forradio interferometry using convex optimization,” in , Pisa, Italy, 2015, pp. 164–168.[52] S. van der Tol, B. D. Jeffs, and A.-J. van der Veen, “Self-calibrationfor the LOFAR radio astronomical array,”
IEEE Transactions on SignalProcessing , vol. 55, no. 9, pp. 4497–4510, 2007.[53] E. Conte, A. De Maio, and G. Ricci, “Recursive estimation of thecovariance matrix of a compound-Gaussian process and its applicationto adaptive CFAR detection,”
IEEE Transactions on Signal Processing ,vol. 50, no. 8, pp. 1908–1915, 2002.[54] F. Pascal, Y. Chitour, J.-P. Ovarlez, P. Forster, and P. Larzabal, “Co-variance structure maximum-likelihood estimates in compound Gaussiannoise: Existence and algorithm analysis,”
IEEE Transactions on SignalProcessing , vol. 56, no. 1, pp. 34–48, 2008.[55] X. Zhang, M. N. El Korso, and M. Pesavento, “Maximum likelihoodand maximum a posteriori direction-of-arrival estimation in the presenceof SIRP noise,” in
International Conference on Acoustics, Speech andSignal Processing , Shanghai, China, 2016.[56] V. Ollier, M. N. El Korso, R. Boyer, P. Larzabal, and M. Pesavento, “Re-laxed concentrated MLE for robust calibration of radio interferometers,”in , Budapest, Hungary,2016.[57] A. Hjorungnes and D. Gesbert, “Complex-valued matrix differentiation:Techniques and key results,”
IEEE Transactions on Signal Processing ,vol. 55, no. 6, pp. 2740–2746, 2007.[58] M. Feder and E. Weinstein, “Parameter estimation of superimposedsignals using the EM algorithm,”
IEEE Transactions on Acoustics,Speech and Signal Processing , vol. 36, no. 4, pp. 477–489, 1988.[59] K. Madsen, H. B. Nielsen, and O. Tingleff, “Methods for non-linear leastsquares problems,” Informatics and Mathematical Modelling, TechnicalUniversity of Denmark, Tech. Rep., 2004.[60] J. Nocedal and S. J. Wright,
Numerical optimization, 2nd ed.
SpringerScience & Business Media, 2006.[61] H. Gavin, “The Levenberg-Marquardt method for nonlinear least squarescurve-fitting problems,” Duke University, 2011.[62] T. W. Anderson,
An introduction to multivariate statistical analysis .Wiley New York, 1958, vol. 2.[63] D. P. Bertsekas,
Nonlinear programming . Athena scientific Belmont,1999.[64] S. Vorobyov, A. B. Gershman, and K. M. Wong, “Maximum likelihooddirection-of-arrival estimation in unknown noise fields using sparsesensor arrays,”
IEEE Trans. Signal Processing , vol. 53, no. 1, pp. 34–43,2005.[65] P. Stoica and R. L. Moses,
Spectral analysis of signals . Pearson PrenticeHall, Upper Saddle River, NJ, 2005.[66] A. Balleri, A. Nehorai, and J. Wang, “Maximum likelihood estimationfor compound-Gaussian clutter with inverse gamma texture,”
IEEETransactions on Aerospace and Electronic Systems , vol. 43, no. 2, pp.775–780, 2007.[67] E. Ollila, D. E. Tyler, V. Koivunen, and H. Poor, “Complex ellipticallysymmetric distributions survey, new results and applications,”
IEEETransactions on Signal Processing , vol. 60, no. 11, pp. 5597–25 625,2012.[68] O. Besson and Y. I. Abramovich, “On the Fisher information matrix formultivariate elliptically contoured distributions,”
IEEE Signal ProcessingLetters , vol. 20, no. 11, pp. 1130–1133, 2013.
69] C. Liu and D. B. Rubin, “ML estimation of the t distribution using EMand its extensions, ECM and ECME,”
Statistica Sinica , vol. 5, no. 1,pp. 19–39, 1995.[70] S. Li, H. Wang, and T. Chai, “A t-distribution based particle filter fortarget tracking,” in
American Control Conference , Minneapolis, MN,2006.
Virginie Ollier was born in Beaumont, France, in1991. She graduated from Ecole Centrale Marseilleand received a Master Research degree in Signaland Image Processing from Ecole Centrale Marseillein 2015. She is currently pursuing a Ph.D. degreein signal processing in University of Paris Saclayand is a member of SATIE laboratory (Syst`emes etApplications des Technologies de l’Information etde l’Energie) and L2S laboratory (Laboratoire desSignaux et Syst`emes). Her research interests includeestimation theory and statistical signal processing,especially for applications in robust array signal processing.
Mohammed Nabil El Korso was born in Oran,Algeria. He received the M.Sc. in Electrical En-gineering from the National Polytechnic School,Algeria in 2007. He obtained the Master Researchdegree in Signal and Image Processing from Paris-Sud XI University, France in 2008. In 2011, heobtained his Ph.D. degree from Paris-Sud XI Univer-sity. From 2011 to 2012, he was a research scientistin the Communication Systems Group at TechnischeUniversit¨at Darmstadt, Germany. He was AssistantProfessor at Ecole Normale Sup´erieure de Cachanfrom 2012 to 2013. Currently, he is Assistant Professor at University of ParisOuest Nanterre la Defense and a member of LEME (EA4416) laboratory.His research interests include statistical signal processing, estimation/detectiontheory with applications to array signal processing.
R´emy Boyer received the M.Sc. and Ph.D degreesfrom the Ecole Nationale Sup´erieure des Telecom-munications (ENST-Paris or T´el´ecom ParisTech) in1999 and 2002, respectively, in statistical signal pro-cessing. From 2002 to 2003, he was a postdoctoralfellow during six months at Sherbrooke University(Canada). From 2003, he is an associate professor atUniversity of Paris-Sud - Laboratory of Signals andSystems (L2S). From 2011 to 2012, R´emy Boyerwas a visiting researcher at the SATIE Laboratory(Ecole Normale Sup´erieure de Cachan) and at theUniversity of Aalborg (Danemark). His teaching activities include basic andadvanced notions of statistical signal processing, in particular in the context ofthe Master2R ATSI (”Automatique et Traitement du Signal et des Images”).R´emy Boyer received an ”Habilitation `a Diriger des Recherches (HDR)” fromUniversity of Paris-Sud in December 2012. His research interests includecompressive sampling, array signal processing, Bayesian performance boundsfor parameter estimation and detection, security in mobile networks as wellas numerical linear and multi-linear algebra.
Pascal Larzabal (M’93) was born in the Basquecountry in the south of France in 1962. He enteredthe Ecole Normale Sup´erieure of Cachan (France), in1985 where he received the Agr´egation in ElectricalEngineering in 1988. He received the PhD in 1992and the ”Habilitation `a Diriger des Recherches” in1998. He is now Professor of Electrical Engineeringat University of ParisSud 11, France. He teacheselectronic, signal processing, control and mathemat-ics. From 1998 to 2003 he was the director of Elec-trical Engineering IUP of the University Paris-Sud.From March 2003 to March 2007 he was at the head of electrical engineeringdepartment in the Institut Universitaire de Technologie of Cachan. SinceJanuary 2007 he is the director of the laboratory SATIE in Paris. His researchconcerns estimation in array processing and spectral analysis for wavefrontidentification, radars, communications, tomography and medical imaging. Hisrecent works concern geographical positioning and radio astronomy.
Marius Pesavento (M’00) received the Dipl.-Ing. and M.Eng. degrees from Ruhr-UniversityBochum, Bochum, Germany, and McMaster Uni-versity, Hamilton, ON, Canada, in 1999 and 2000,respectively, and the Dr.-Ing. degree in electricalengineering from Ruhr-University Bochum in 2005.Between 2005 and 2009, he held research positionsin two start-up companies. In 2010, he became anAssistant Professor for Robust Signal Processingand a Full Professor for Communication Systems in2013, at the Department of Electrical Engineeringand Information Technology, Technical University Darmstadt, Darmstadt,Germany. His research interests include robust signal processing and adap-tive beamforming, high-resolution sensor array processing, multiantenna andmultiuser communication systems, distributed, sparse, and mixed-integeroptimization techniques for signal processing and communications, statisticalsignal processing, spectral analysis, and parameter estimation. He has receivedthe 2003 ITG/VDE Best Paper Award, the 2005 Young Author Best PaperAward of the IEEE Transactions on Signal Processing, and the 2010 BestPaper Award of the CrownCOM conference. He is a Member of the EditorialBoard of the EURASIP Signal Processing Journal, and served as an AssociateEditor for the IEEE Transactions on Signal Processing in 2012-2016. He isa Member of the Sensor Array and Multichannel Technical Committee ofthe IEEE Signal Processing Society, and the Special Area Teams “SignalProcessing for Communications and Networking” and “Signal Processing forMultisensor Systems” of the EURASIP.(M’00) received the Dipl.-Ing. and M.Eng. degrees from Ruhr-UniversityBochum, Bochum, Germany, and McMaster Uni-versity, Hamilton, ON, Canada, in 1999 and 2000,respectively, and the Dr.-Ing. degree in electricalengineering from Ruhr-University Bochum in 2005.Between 2005 and 2009, he held research positionsin two start-up companies. In 2010, he became anAssistant Professor for Robust Signal Processingand a Full Professor for Communication Systems in2013, at the Department of Electrical Engineeringand Information Technology, Technical University Darmstadt, Darmstadt,Germany. His research interests include robust signal processing and adap-tive beamforming, high-resolution sensor array processing, multiantenna andmultiuser communication systems, distributed, sparse, and mixed-integeroptimization techniques for signal processing and communications, statisticalsignal processing, spectral analysis, and parameter estimation. He has receivedthe 2003 ITG/VDE Best Paper Award, the 2005 Young Author Best PaperAward of the IEEE Transactions on Signal Processing, and the 2010 BestPaper Award of the CrownCOM conference. He is a Member of the EditorialBoard of the EURASIP Signal Processing Journal, and served as an AssociateEditor for the IEEE Transactions on Signal Processing in 2012-2016. He isa Member of the Sensor Array and Multichannel Technical Committee ofthe IEEE Signal Processing Society, and the Special Area Teams “SignalProcessing for Communications and Networking” and “Signal Processing forMultisensor Systems” of the EURASIP.