Calibration and Direction-of-Arrival (DoA) estimation is a major issue in array processing [2, 3]. The latter has been studied in several applications, e.g., radar, sonar, satellite, wireless communication and radio interferometric systems [4, 5], where we commonly use largely distributed sensors elements aiming to achieve high resolution. In all these sensor network applications, calibration is required as some parameters are not exactly known due to imperfect instrumentation or propagation conditions [6]. Let us note that calibration algorithms are distinguished by the presence [7] or absence [8] of one or more cooperative sources, named calibrator sources. Indeed, prior source information can be available [6] and consists mainly in the true/nominal directions and powers of calibrator sources (i.e., without any perturbation effects or antenna imperfections). Furthermore, most calibration algorithms are based on the least squares approach, with a sequential procedure updating each parameter alternatively [4]. The least squares estimator is indeed equivalent to the Maximum Likelihood (ML) method under a (unrealistic) Gaussian noise model.
The aim of the proposed methodology here is to estimate successively the unknown sensor gains and phase errors, along with the calibrator and noise parameters, through minimization of a proper weighting cost function. In this work, uncertainties are estimated from the array covariance matrix, since dealing directly with time series data and operating on the signal domain quickly becomes computationally unfeasible for a large number of samples [9]. The scenario under study is general but could be adapted to any practical application as in the radio astronomy context, where the number of parameters to estimate is tremendous and frequency bands are wide.
In the multi-frequency scenario, a suboptimal way to perform calibration is to consider one wavalength bin at a time, with only one centralized processor, which has access to data in the whole available range of wavelengths. In this work, we study an accelerated version based on the scalable form of the Alternating Direction Method of Multipliers (ADMM) [10, 11] with a specific network topology: there is no fusion center and agents exchange information only among themselves. The goal being to reduce the complexity in operation flow and signaling exchanging [12, 13, 14, 15, 16, 17]. For estimation of the directional gains, the compressive sensing framework, especially the sparse representation method, is well-adapted and has already been applied for source localization in fully and partially calibrated arrays [18, 19, 20, 21].
The notation used through this paper is the following: denote, respectively, the complex conjugate, transpose, Hermitian operator, element-wise raising to
, real part and the n-th element of a vector. The expectation operator is
respectively, the Kronecker, the Khatri-Rao and the Hadamard product. The operator diag(.) converts a vector into a diagonal matrix, blkdiag(.) is the block-diagonal operator, whereas vecdiag(.) produces a vector from the main diagonal of a matrix and vec(.) stacks the columns of a matrix on top of one another. The operators
refer to the
and Frobenius norms, respectively. Finally,
identity matrix and
refers to the cardinality of a set.
Let us consider Q emitting signal sources and P sensor elements in the array. Each source direction is defined by a 2-dimensional vector
, s.t., all nominal/true known directions, without any
disturbances, are stacked in Propagation conditions induce wavelength dependent distortions, leading to apparent source directions
different from the true ones. Under the narrowband assumption, the array response matrix reads
includes the known Cartesian coordinates describing each sensor location in the array, s.t., for
. Therefore, the
narrowband signals measured by all antennas is written as follows, for the n-th time sample and wavelength
where the undirectional antenna gains are collected in the complex diagonal matrix and the directional gain responses, assumed identical for all antennas, are modeled by the diagonal matrix
Finally,
are the i.i.d. calibrator source signal and additive Gaussian thermal noise vectors with their corresponding diagonal covariance matrices
, respectively. From (1), we deduce the following covariance matrix
where In this context, the calibration problem consists in estimating the parameter vector of interest
the total number of available wavelengths and
. To this end, we exploit sample co- variance matrices
, defined as
for wavelength
In estimation theory, the ML estimator is well-known for its statistical efficiency but not always easy to implement in practice. The Weighting Least Squares approach is an appropriate alternative as it is asymptotically equivalent to the ML for a large number of samples N. Therefore, we wish to minimize the following local cost function, associated to wavelength
where . Most sources are assumed buried beneath the noise and antennas are identical in the array with negligible mutual coupling. The aim of the designed calibration algorithm is to minimize the global cost function
in a parallel and step-wise approach, with
the total set of available wavelengths. Usually, minimization is conducted w.r.t. one specific parameter while fixing the others in
[22].Here, our approach is different: we propose an accelerated version where estimation is performed directly w.r.t. the consensus (hidden) variables, as described in Algorithm 1 and detailed in the following.
To achieve multi-frequency calibration in the sensor array, coherence is imposed along wavelength subbands for both directional and undirectional gains, by imposing available constraints or enforcing smooth variation. The choice of the basis functions is motivated by the application under analysis and can be adapted accordingly.
3.1 Coherence model for the undirectional antenna gains
To impose coherence along subbands, we introduce a set of smooth wavelength dependent basis functions and express the gains as linear combinations. Let us define , the consensus vector for the p-th sensor with unknown linear coefficients. Therefore, for
stands for the polynomial terms, describing the variation of the undirectional gains w.r.t. wavelength. For instance, we
can consider the typical basis function is the studied frequency of interest with c the speed of light and
reference frequency [22, 23]. By stacking all vectors
, we obtain the global consensus vector
leading to
with
3.2 Coherence model for the directional gains
Similarly as for the undirectional gains, the coherence model is defined as follows: let us consider , such that for
in which is the vector of hidden variables for the q-th calibrator source, associated to directional gains
is the corresponding basis vector. As in section 3.1, all
are stacked in
, finally leading to
with . We assume identical behavior for all sources but the process can be straightforwardly adapted to different behavior. In [22], the directional gains in
were assumed inversely proportional to
the algorithm can be adjusted to any general existing models.
3.3 Distributed network with a fusion center
Dealing with large data volumes delivered by advanced sensor array systems requires computationally efficient calibration algorithms, with a huge number of unknowns to solve. To improve both computational cost and estimation accuracy, distributed calibration has been proposed by exploiting data parallelism across frequency. Contrary to a centralized hardware architecture which processes all frequency bands at a single location and is therefore computationally challenging, distributed optimization introduces more than one compute agents and analyzes the data simultaneously across smaller frequency intervals [10]. By distributing the total computations across the network, we gain a significant reduction in operational and energy cost and each agent receives information indirectly across the whole frequency range, thus improving the calibration accuracy. To handle this, let us consider Z computational agents disposed on a network. Each agent has access to some wavelengths . The corresponding unknown pa- rameters in p are estimated locally and consensus is enforced among agents by imposing constraints in (4) and (6).
To start with, let us focus on estimation of the undirectionnal sensor gains in section 3.1. We define as the local copy of the common optimization variable
-th agent and we note
the set of all
in the network. Calibration is reformulated as the following constrained problem
where is the cost function for the z-th agent, i.e., for
depends on the local variable
and is associated to data
To solve this problem, we use the augmented Lagrangian, given by [24]
where Lagrange multipliers and
is the regularization term. We resort to the consensus ADMM in the scaled form by introducing the scaled dual variable
[10]. The three updates of the iterative
algorithm are therefore given by
where t is the iteration counter. Minimization (9) leads to the following average, computed at the fusion center and sent to all agents in the network,
from which the undirectional gains can be directly deduced with (4). The local minimization step in (8) is the computationally most expensive one. To this end, we adopt an iterative approach and notice that the problem is separable w.r.t. each , i.e., w.r.t. each agent. Let us assume
as two independent variables [25]. We then minimize
w.r.t.
, considering
as fixed and neglecting the diagonal elements in the cost function. In this case, the local cost function becomes separable w.r.t. the sub-vectors of
local consensus vector for the p-th sensor at the z-th agent. The following decompositions w.r.t. the sensor elements are also possible
and ρ
corresponds to the cost function for the p-th row of
, which only depends on
since the remaining parameters are considered as fixed in this step. Let us define the operator
converts to a vector the p-th row of a matrix and removes the p-th element of this selected vector. We also introduce the quantity
(reference source model) and the following vectors
in which In addition, let us consider the
in (12) as
and finally obtain the following estimate
3.4 Distributed network with no fusion center
We consider a specific formulation of the ADMM where every node in the network performs calibration locally and consensus is only reached with clearly identified neighbours without fusion center [12]. We note the index set that corresponds to the neighbours of the z-th agent. The considered network architecture is exposed in Figure 1 where for example,
define the quantity
as the copy available at the z-th agent, transferred to the y-th agent. In such context, the minimization problem becomes
where the auxiliary variables impose consensus contraints on two neighboring agents and are meant to be local copies of
. The decentralized strategy enables to cooperatively minimize a sum of local objective functions, the final aim being to converge to a common value, with fast convergence speed and good estimation performance [26]. To obtain a more compact form of the problem in (15), we define
leading to
with . As in section 3.3, the scaled version of the ADMM leads to
and through decomposition of the problem in (17) w.r.t. sensor dependence, we obtain
with . The selected variables
are obtained from
via an appropriate selection matrix. After considering the projection onto B and denoting the messages passed between the agent as
we solve (18) thanks to
The steps of the proposed distributed method for calibration of sensor gains are exposed in Algorithm 1.2.
3.5 Estimation of directional gains
In this section, we describe the part of the algorithm dedicated to the estimation of DoA and directional gains
, for fixed sensor gains, with a sparse and distributed implementation. Assuming a sparse observed scene, we define dictionaries of steering matrices for
notes the total number of directions on the grid. The sparse vectors in contain the corresponding squared direc-
tion dependent gains. The covariance model is rewritten as
estimation and satisfy both sparsity and positivity requirements, we use the Distributed Iterative Hard Thresholding (IHT) [27, 28]. But contrary to [22],
the following hard-thresholding operator
ered to provide access to the DoA of the q-th source, and a first estimate of the directional gain . The quantity
refers to the q-th column of a ma- trix, the expression
discards the elements corresponding to the diagonal of
and the hard thresholding operator
-largest components and sets the remaining entries equal to zero. Finally, thanks to (5) and dealing with the consensus variables as in section 3.4, the minimization problem becomes
where we benefit from the previous hard-thresholding estimate to define
. As previously, we impose consensus between neighbours thanks to some auxiliary variables but due to lack of space, we only present here the resulting local update for
where , we obtain an estimate of
and process the next source, as shown in Algorithm 1.3.
In order to evaluate the method, we consider realistic simulations for the radio astronomy context where the new generation of phased array systems such as the Low Frequency Array (LOFAR) and the Square Kilometre Array (SKA) requires the development of new advanced signal processing techniques for calibration purpose [4, 29]. Indeed, lack of calibration leads to dramatic effects and distortions in the reconstructed images. We consider P = 60 antennas spread over a five-armed spiral [30, 31], which corresponds to the LOFAR’s Initial Test Station. Let us assume a sky model with Q = 3 strong calibrator sources and weak unknown sources in the background. The reference frequency
and we consider frequencies ranging from 29.6 MHz to 30.4 MHz, with Z = 3 agents in the network and
. The polynomial orders are chosen as
The consensus variables
are initialized as zeros and the squared directional gains are generated thanks to power law functions
4.1 Influence of the number of frequency channels
First of all, we investigate the statistical performance of the proposed distributed algorithm as a function of the number of samples N or the Signal-to-Noise Ratio (SNR). The SNR is defined as the ratio between the sum of apparent powers for all Q sources and the noise power. Results are averaged for 100 Monte-Carlo runs. In Figure 2, we plot the three following cases: F = 3 and each agent handles one frequency, i.e., (green curve),
(blue curve) and
(red curve). In Figure 2 (a), we plot the Root Mean Square Error (RMSE) as a function of N for the undirectional gains
, defined as
for fixed SNR =
dB. A similar figure is presented in Figure 2 (b), for the source directions
, as a function of the SNR and fixed
illustrate the performance by comparing with the mono-calibration scenario where each agent handles one single frequency, independently. We notice that mono-calibration is clearly improved, by using a distributed procedure where the whole information is flowing through the entire network.
4.2 Influence of the network architecture
We aim to show the advantages of the proposed distributed network with no fusion center and only exchange of local information between neighboring agents, in terms of complexity. With similar number of iterations in all loops of the algorithm, different estimation performance are attained in Figure 2 (a) while similar RMSE is reachable in Figure 2 (b) but with an additional computational cost if there is a fusion center (an increase of at least a factor 5 in computing time).
4.3 Convergence analysis
We illustrate the convergence behavior of the proposed algorithm by analyzing the following residuals as function of the iteration number. Depending on the iteration in Algorithm 1, we plot the primal residual as a function of the iteration number of Algorithm 1.2, defined as
Likewise, we also study the different estimates between agents through
Similar statistical behavior can be ob- tained for corresponding residuals in Algorithm 1.3
In this work, we proposed an iterative algorithm for parallel calibration, applied in a general context of sensor array processing: complex electronic gains are imprecisely known and propagation disturbances lead to deviations in the source locations. In order to reduce the communication overhead, the spe-cific variation of parameters across wavelength is exploited in a distributed network with no fusion center and local exchange of information between adjacent connected nodes. The two main steps of the algorithm are based on the scalable form of the ADMM and distributed IHT procedures. We highlighted the effectiveness and time efficiency of the proposed method using simulated data, even in the presence of non-calibrator sources at unknown directions.
Figure 1: Example of distributed network with no fusion center.
[1] M. Brossard, V. Ollier, M. N. El Korso, R. Boyer, P. Larzabal, Fast decentralized multi-frequency calibration for doa estimation with distributed sensors, in: DAT 2020, Proceeding Conf. Algiers, Algeria, 2020.
[2] P. Stoica, K. C. Sharman, Maximum likelihood methods for direction-of- arrival estimation, IEEE Transactions on Acoustics, Speech and Signal Processing 38 (7) (1990) 1132–1143.
[3] S. Vorobyov, A. B. Gershman, K. M. Wong, Maximum likelihood direction-of-arrival estimation in unknown noise fields using sparse sensor arrays, IEEE Transactions on Signal Processing 53 (1) (2005) 34–43.
[4] A.-J. van der Veen, S. J. Wijnholds, Signal processing tools for radio astronomy, in: Handbook of Signal Processing Systems, Springer, 2013, pp. 421–463.
[5] L. C. Godara, Application of antenna arrays to mobile communications, part II: Beam-forming and direction-of-arrival considerations, Proceedings of the IEEE 85 (8) (1997) 1195–1245.
[6] B. C. Ng, C. M. S. See, Sensor-array calibration using a maximum- likelihood approach, IEEE Transactions on Antennas and Propagation 44 (6) (1996) 827–835.
[7] B. C. Ng, A. Nehorai, Active array sensor localization, Elsevier Signal processing 44 (3) (1995) 309–327.
[8] A. J. Weiss, B. Friedlander, Array shape calibration using sources in unknown locations-a maximum likelihood approach, IEEE Transactions on Acoustics, Speech, and Signal Processing 37 (12) (1989) 1958–1966.
[9] S. J. Wijnholds, S. Chiarucci, Blind calibration of phased arrays us- ing sparsity constraints on the signal model, in: 24th European Signal Processing Conference (EUSIPCO), 2016, pp. 270–274.
[10] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, Distributed opti- mization and statistical learning via the alternating direction method of multipliers, Foundations and Trends in Machine Learning 3 (1) (2011) 1–122.
[11] S. Kazemi, S. Yatawatta, S. Zaroubi, Clustered calibration: an im- provement to radio interferometric direction-dependent self-calibration, Monthly Notices of the Royal Astronomical Society 430 (2) (2013) 1457– 1472.
[12] T. Erseghe, A distributed and scalable processing method based upon ADMM, IEEE Signal Processing Letters 19 (9) (2012) 563–566.
[13] W. Shi, Q. Ling, K. Yuan, G. Wu, W. Yin, On the linear convergence of the ADMM in decentralized consensus optimization., IEEE Transactions on Signal Processing 62 (7) (2014) 1750–1761.
[14] J. F. C. Mota, J. M. F. Xavier, P. M. Q. Aguiar, M. Puschel, D-ADMM: a communication-efficient distributed algorithm for separable optimization, IEEE Transactions on Signal Processing 61 (10) (2013) 2718–2723.
[15] V. Ollier, M. N. El Korso, A. Ferrari, R. Boyer, P. Larzabal, Robust distributed calibration of radio interferometers with direction dependent distortions, accepted for publication by Elsevier Signal Processing.
[16] S. Yatawatta, Fine tuning consensus optimization for distributed radio interferometric calibration, in: 24th European Signal Processing Conference (EUSIPCO), Budapest, Hungary, 2016.
Figure 2: (a) RMSE on the undirectional gains as function of the number of samples N, (b) RMSE on the apparent source directions as function of the SNR.
[17] S. Yatawatta, S. Kazemi, S. Zaroubi, GPU accelerated nonlinear op- timization in radio interferometric calibration, in: Innovative Parallel Computing (InPar), IEEE, San Jose, CA, 2012, pp. 1–6.
[18] D. Malioutov, M. Cetin, A. S. Willsky, A sparse signal reconstruction perspective for source localization with sensor arrays, IEEE Transac-
Figure 3: (a) Statistical comparison between different network topologies for same computational cost, (b) Statistical comparison between different network topologies for different computational cost.
tions on Signal Processing 53 (8) (2005) 3010–3022.
[19] C. Steffens, P. Parvazi, M. Pesavento, Direction finding and array cal- ibration based on sparse reconstruction in partly calibrated arrays, in: 8th Sensor Array and Multichannel Signal Processing Workshop (IEEE SAM), 2014, pp. 21–24.
Figure 4: (a) Primal residual and (b) estimates difference
of the local consensus variables among agents as function of the iteration t in Algorithm 2, for different values of the iteration i in Algorithm 1.
[20] M. Haardt, M. Pesavento, F. Roemer, M. N. El Korso, Subspace meth- ods and exploitation of special array structures, in: Electronic Reference in Signal Processing:Array and Statistical Signal Processing (M. Virberg, ed.), Vol. 3, Academic Press Library in Signal Processing, Elsevier Ltd., 2014.
[21] V. Ollier, M. N. El Korso, R. Boyer, P. Larzabal, M. Pesavento, Joint ML calibration and DOA estimation with separated arrays, in: International Conference on Acoustics, Speech and Signal Processing (IEEE ICASSP), Shanghai, China, 2016.
[22] M. Brossard, M. N. El Korso, M. Pesavento, R. Boyer, P. Larzabal, S. J. Wijnholds, Parallel multi-wavelength calibration algorithm for radio astronomical arrays, Elsevier Signal Processing Journal 145 (2018) 258–271.
[23] S. Yatawatta, Distributed radio interferometric calibration, Monthly Notices of the Royal Astronomical Society 449 (4) (2015) 4506–4514.
[24] L. Li, X. Wang, G. Wang, Alternating Direction Method of Multipliers for separable convex optimization of real functions in complex variables, Mathematical Problems in Engineering 2015 (Article ID 104531) (2015) 1–14.
[25] S. Salvini, S. J. Wijnholds, Fast gain calibration in radio astronomy using alternating direction implicit methods: Analysis and applications, Astronomy & Astrophysics 571 (A97).
[26] T. Erseghe, D. Zennaro, E. Dall’Anese, L. Vangelista, Fast consensus by the alternating direction multipliers method, IEEE Transactions on Signal Processing 59 (11) (2011) 5523–5537.
[27] T. Blumensath, M. E. Davies, Normalized iterative hard thresholding: Guaranteed stability and performance, IEEE Journal of selected topics in signal processing 4 (2) (2010) 298–309.
[28] S. Patterson, Y. C. Eldar, I. Keidar, Distributed compressed sensing for static and time-varying networks, IEEE Transactions on Signal Processing 62 (19) (2014) 4931–4946.
[29] V. Ollier, M. N. El Korso, R. Boyer, P. Larzabal, M. Pesavento, Relaxed concentrated MLE for robust calibration of radio interferometers, in: 24th European Signal Processing Conference (EUSIPCO), Budapest, Hungary, 2016, pp. 280–284.
[30] S. J. Wijnholds, J. D. Bregman, A. J. Boonstra, Sky noise limited snap- shot imaging in the presence of RFI with Lofar’s Initial Test Station, Experimental Astronomy 17 (1-3) (2004) 35–42.
[31] V. Ollier, M. N. El Korso, R. Boyer, P. Larzabal, M. Pesavento, Robust calibration of radio interferometers in non-Gaussian environment, IEEE Transactions on Signal Processing 65 (21) (2017) 5649–5660.
[32] E. Ghadimi, A. Teixeira, I. Shames, M. Johansson, Optimal parameter selection for the alternating direction method of multipliers (ADMM): quadratic problems, IEEE Transactions on Automatic Control 60 (3) (2015) 644–658.