Noncoherent Sensor Fusion via
entropy regularized
optimal mass transport
Abstract
This work presents a method for information fusion in source localization applications. The method utilizes the concept of optimal mass transport in order to construct estimates of the spatial spectrum using a convex barycenter formulation. We introduce an entropy regularization term to the convex objective, which allows for lowcomplexity iterations of the solution algorithm and thus makes the proposed method applicable also to higherdimensional problems. We illustrate the proposed method’s inherent robustness to misalignment and miscalibration of the sensor arrays using numerical examples of localization in two dimensions.
Noncoherent Sensor Fusion via
entropy regularized optimal mass transport
Filip Elvander, Isabel Haasler, Andreas Jakobsson, and Johan Karlsson^{†}^{†}thanks: This work was supported in part by the Swedish Research Council, Carl Trygger’s foundation, and the Royal Physiographic Society in Lund. 
Div. of Mathematical Statistics, Lund University, Sweden 
Dept. of Mathematics, KTH Royal Institute of Technology, Sweden 
emails: {filipelv, aj}@maths.lth.se, {haasler, johan.karlsson}@math.kth.se 
Index Terms— Optimal mass transport, Entropy regularization, Target localization, Sensor fusion, Noncoherent processing.
1 Introduction
Source localization using an array of sensors is a well established and still highly active field of research in signal processing, with applications in areas such as radar, sonar, and audio processing [1, 2, 3, 4]. The localization problem is often formulated as a spatial spectral estimation problem, where signal sources are identified with peaks of power in the spectral representation. In contrast to the closely related problem of direction of arrival (DoA) estimation, where one aims to find the directions of incoming planar waves from far field sources, the problem of localization also includes the estimation of the range (i.e., distance to the sources) and thus needs to utilize that the incoming wavefronts are nonplanar. Classical methods, originally used mainly for DoA estimation, such as beamforming techniques [5] and subspace methods like MUSIC [6] and ESPRIT [7], may be readily extended to the localization problem, with more recent efforts including works exploiting spectral sparsity [8, 9, 10]. Many of these methods rely on the assumption of wellcalibrated arrays, i.e., known sensor positions, in order to produce reliable estimates. However, this assumption may in some applications be violated due to, e.g., calibration errors, perturbations of the sensor localizations, or synchronization errors. Therefore, considerable effort has been directed towards formulating estimators robust to such calibration errors (see, e.g., [11, 12, 13, 14, 15]). In some of the approaches to address this problems, techniques have been developed that divide the available sensors into subarrays, each assumed to be perfectly calibrated, although relative errors between the different subarrays are allowed for [16]. The DoA estimation and/or localization is then performed by minimizing measures of discrepancy between the measurements, i.e., the estimated joint covariance matrix of the complete set of sensors, and the assumed array geometry [15].
In this work, we instead consider localization by fusing information from separate sensor arrays, where we only assume to have access to measurements in the form of covariance matrices for the individual arrays, socalled noncoherent processing (see, e.g., [17, 18, 19, 20]), i.e., we assume no information about crosscorrelation between different arrays. Such a scenario may arise in distributed sensor networks, wherein sensor nodes transmit information to a central processing unit that performs the estimation [21, 22]. In such settings, it is common to have restrictions on the communication bandwidth, which requires transmitting only aggregated data [21], e.g., in the form of covariance matrices. Also, as the separate sensor arrays may not be jointly synchronized, crosscorrelations cannot be estimated from the raw data.
Extending on the work presented in [23], which proposes a robust method for DoA sensor fusion using optimal mass transport (OMT), we here generalize this approach to the source localization problem. One of the main ideas is to use OMT as a notion of distance between spectra (see, e.g., [24, 25, 26]), and in particular, the proposed method performs sensor fusion by computing the corresponding barycenter solutions based on incomplete measurements. In addition, to make the scheme computationally efficient, we introduce an entropy regularization term in the problem formulation and show how to solve the resulting optimization problem iteratively, similar to the popular Sinkhorn iteration scheme [27] (see also [28]). The robustness of the proposed method is demonstrated in numerical examples.
2 Background
2.1 Near field localization
Consider an array consisting of sensors located at , for , and assume that pointlike sources, located in a space , emit spherical waves impinging on the sensors. Then, the spatial distribution of the energy of the source signals may be representated by a nonnegative function, or measure, on , referred to as the power spectrum. We will denote the space of such functions by . Letting denote the wavelength of the impinging waves, the covariance matrix of the sensor measurements is , where the linear operator is given by (see, e.g., [29])
(1) 
where is the array manifold vector [30], i.e.,
(2) 
with denoting the Lebesgue measure. In other words, a positive semidefinite matrix is a valid covariance for a sensor array if and only if , for some . From this, it may be noted that performing source localization based on sensor measurements corresponds to inferring a power spectrum from observations of the array covariance matrix .
We will here consider a scenario wherein several sensor arrays measure the superposition of the impinging waves, i.e., when covariance matrices , for , are available, corresponding to separate sensor arrays. For the case when all arrays are calibrated, i.e., when the position of the sensors are known perfectly, source localization thus corresponds to finding a single power spectrum that is consistent with the observed matrices, i.e., , for , where is the linear operator corresponding to array . However, for the case of miscalibration, the (assumed) operators contain errors, which may result in there being no single consistent with all covariance matrices. We will here adress this problem by utilizing OMT in order to induce robustness to the sensor fusion problem.
2.2 Optimal mass transport
Let and be two nonnegative mass distributions on the space . A transport plan for such a pair is a nonnegative distribution on , with , where describes the amount of mass transported from location to location . By defining a cost function describing the cost of moving mass, we may define the total cost associated with a transport plan as
(3) 
The OMT problem is then to find the feasible transport plan from to incurring minimal transportation cost (see, e.g., [31]), i.e., to solve
(4)  
subject to  
where the constraints ensure that the transport plan is valid, i.e., matches the marginals and . The suitability of utilizing OMT costs as measures of distance has previously been considered for imposing metric structure on the space of power spectra [24], as well as to induce distances and smooth interpolants for covariance matrices [25]. In this work, the optimal transport cost will be used for interpolating between mass distributions and . It is worth noting that the OMT cost quantifies location errors based on the distances on the underlying space, in this case , thereby allowing for offsets due to, e.g., calibration errors. This is in contrast to standard metrics (e.g., , ) which only compare functions point by point without taking the geometry of the underlying space, i.e., relative location of the points, into account.
3 Proposed method
Consider the localization scenario with the source signals impinging on a set of arrays, with array manifold vectors , generating the covariance matrices , . Given the array geometries and covariance matrix estimates of each individual array, we aim to estimate the spatial spectrum . That is, we assume that each sensor array is calibrated within itself, but there might be global calibration errors resulting in inaccurate knowledge of each sensor array’s position and orientation. For the far field case, i.e., DoA estimation, an OMT formulation for fusing information from different sensor arrays has been proposed in [23]. There, the estimate was formed as the spectrum closest, in the OMT sense, to a set of spectra, each consistent with a corresponding covariance matrix. Correspondingly, in the localization setting the spatial spectrum may be found as the solution to the convex barycenter problem
(5)  
subject to 
where denotes the linear operator (1) of array . In order to find practical solutions to the barycenter problem in (5), the space may be discretized as to be represented by points . Accordingly, the marginal spatial spectra are represented by the vectors , and the cost function and transport plans are represented by matrices and , respectively, where , and denotes the mass transported from to . Although such a discretization yields finitedimensional linear programs, finding an optimal transport plan can be computationally demanding due to the large number of variables.
For discretizations of the standard OMT problem (4), this was in [27] addressed by regularizing the problem by introducing an entropy term to the objective function, yielding the strictly convex program
(6)  
subject to 
where is an entropy term, a small constant, and denotes a vector of ones. In [27], it was shown that this problem may be solved in a computationally efficient way, using socalled Sinkhorn iterations.
3.1 Efficient sensor fusion
In order to arrive at a computationally efficient scheme for computing spectral barycenters, we propose to regularize the problem in (5) according to
(7)  
subject to 
Here, , for , are the vectorizations of the covariances matrices, and are matrix representations of the discretized operators . Also, the vectors allow for noisy covariance matrix estimates, penalized by the standard norm, as controlled by the positive penalty parameter . The addition of the entropy term, , renders the problem strictly convex and allows for efficiently computing the barycenter , using the observation in [28], that Sinkhorn iterations are equivalent to a block coordinate ascent in a Lagrange dual of the entropy regularized transport problem. This then allows for increasing the dimensionality of the sensor fusion problem, i.e., making it feasible to include also range estimation or considering higher spatial dimensions. The extended Sinkhorn iterations, constituting a globally convergent solution algorithm for (7), are detailed in the following proposition.
Proposition 1
A block coordinate ascent scheme for the dual of the barycenter problem (7) is to iteratively perform

For , let be the solution to
(8) where

Update
where
Here, , and , , and denotes the elementwise exponential function, elementwise division, and elementwise product, respectively. The barycenter solution is then given by , where is the limit point of the iterations. Note that we here convert all complex quantities to an equivalent realvalued representation such that the vectorized covariances, , are mapped according to
where and denote the real and imaginary parts, respectively, and where and the discretized operators are converted as to be consistent with this.
A proof of the proposition may be found in Appendix A. It may be noted that for the noisefree case with complete observations of the underlying spectra, i.e., and , as well as , one obtains the standard Sinkhorn iterations. The roots of (8) are determined by the means of a Newton method, using the Jacobian
The complete estimation method for computing the spectral barycenter is summarized in Algorithm 1.
4 Numerical results
In this section, we illustrate the proposed method’s behavior in sensor fusion scenarios for localization in 2D, demonstrating the robustness to array alignment errors that the OMT criterium induces. Specifically, we consider two uncorrelated sources impinging on two asynchronous sensor arrays; one ellipsoidal shaped array consisting of 8 sensors, and one linear array consisting of 7 sensors. To simulate misalignment between the two arrays, we introduce an unknown rotation to the ellipsoidal array. The scenario, as shown in Figure 1, is then used in a simulation study. Varying the misalignment, i.e., the rotation angle, between 0 and 10 degrees, we for each considered angle perform 100 Monte Carlo simulations, were the locations of the signal sources are randomized uniformly on the square . We let the source signals be uncorrelated circularly symmetric Gaussian white noises and add an uncorrelated Gaussian noise to the sensors. The source signals have variance 100 and the variance of the sensor noise is 1. The wavelength of the impinging waves is twice that of the smallest sensor spacing in the linear array. Estimating the covariance matrices of the separate sensor arrays using the sample covariance matrix and 500 signal snapshots, we apply the proposed estimator in (7) to the covariance matrix estimates and record the distance from the location estimates to the actual target positions. Throughout, we use the parameters and , and the cost function , for grid points .
As comparison, we also obtain estimates using the noncoherent MUSIC and MVDR estimators, as described in [18], as well as the leastsquares (LS) estimator from [32], and the noncoherent SPICE estimator from [19]. Figures 1 and 2 show the estimates obtained in one simulation using the proposed and MUSIC estimators, respectively. The misalignment angle is 6.7 degrees. As can be seen, for the proposed estimator, the misalignment only causes a small perturbation on the estimated source locations, demonstrating the ability of the OMT formulation to perform smooth interpolation on the underlying domain. In contrast, the MUSIC (pseudo) spectral estimate contains spurious peaks, as well as larger deviations from the true source locations. The robustness of the proposed estimator is further demonstrated in Figure 3, displaying the results from the simulations. From this, it may be noted that the proposed estimator deviates slightly more from the ground truth for the nonperturbed case as compared to the other methods. However, the robustness to misalignment is significantly higher; the rate of error increase is considerably lower than for the other methods.
Appendix A Proof of Proposition
Proof 1
In the following, we denote the set by , and similarly for , and . The Lagrangian of (7) with the dual variables and is
(9) 
Minimizing the Lagrangian with respect to gives the expression for the mass transport matrix solutions,
(10) 
with , and , for . Further, minimizing (9) with respect to the variables and , one gets a dual problem to (7), formulated as
(11)  
subject to 
where
A block coordinate ascent method is to maximize (11) iteratively with respect to and . Maximization with respect to requires
which gives the first step (8). To maximize (11) with respect to , consider the corresponding Lagrangian with dual variable , given by
Maximization with respect to requires
Further, due to the constraint in (11), it follows
This completes the second step of the method. Finally, with (10), we get the expression for the barycenter marginal
References
 [1] H. Krim and M. Viberg, “Two Decades of Array Signal Processing Research,” IEEE Signal Process. Mag., pp. 67–94, July 1996.
 [2] J. C. Chen, Kung Yao, and R. E. Hudson, “Source localization and beamforming,” IEEE Signal Process. Mag., vol. 19, no. 2, pp. 30–39, Aug. 2002.
 [3] P. Stoica and R. Moses, Spectral Analysis of Signals, Prentice Hall, Upper Saddle River, N.J., 2005.
 [4] N. Dey and A. S. Ashour, Direction of Arrival Estimation and Localization of MultiSpeech Sources, Springer, 2018.
 [5] J. Capon, “High Resolution Frequency Wave Number Spectrum Analysis,” Proc. IEEE, vol. 57, pp. 1408–1418, 1969.
 [6] R. Schmidt, “Multiple emitter location and signal parameter estimation,” in Proceedings of RADC Spectrum Estimation Workshop, 1979, pp. 243–258.
 [7] R. Roy and T. Kailath, “ESPRIT – Estimation of Signal Parameters via Rotational Invariance Techniques,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, no. 7, pp. 984–995, July 1989.
 [8] P. Stoica, P. Babu, and J. Li, “SPICE : a novel covariancebased sparse estimation method for array processing,” IEEE Trans. Signal Process., vol. 59, no. 2, pp. 629 –638, Feb. 2011.
 [9] S. I. Adalbjörnsson, T. Kronvall, S. Burgess, K. Åström, and A. Jakobsson, “Sparse Localization of Harmonic Audio Sources,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 24, no. 1, pp. 117–129, January 2016.
 [10] Z. Yang and L. Xie, “Enhancing Sparsity and Resolution via Reweighted Atomic Norm Minimization,” IEEE Trans. Signal Process., vol. 64, no. 4, pp. 995–1006, Feb 2016.
 [11] J. Yang and A. Swindlehurst, “The Effects of Array Calibration Errors on DFBased Signal Copy Performance,” IEEE Trans. Signal Process., vol. 43, no. 11, pp. 2724–2732, November 1995.
 [12] S. D. Somasundaram and A. Jakobsson, “Degradation of Covariance ReconstructionBased Robust Adaptive Beamformers,” in Sensor Signal Processing for Defense, Edinburgh, Sept. 89 2014.
 [13] M. Rübsamen and M. Pesavento, “Maximally Robust Capon Beamformer,” IEEE Trans. Signal Process., vol. 61, no. 8, pp. 2030–2041, 2013.
 [14] R. G. Lorenz and S. Boyd, “Robust Minimum Variance Beamforming,” IEEE Trans. Signal Process., vol. 53, no. 5, pp. 1684–1696, May 2005.
 [15] M. Pesavento, A. B. Gershman, and K. M. Wong, “Direction finding in partly calibrated sensor arrays composed of multiple subarrays,” IEEE Trans. Signal Process., vol. 50, no. 9, pp. 2103–2115, Sep. 2002.
 [16] C. M. See and A. B. Gershman, “Directionofarrival estimation in partly calibrated subarraybased sensor arrays,” IEEE Trans. Signal Process., 2004.
 [17] M. Wax and T. Kailath, “Detection of Signals by Information Theoretic Criteria,” IEEE Trans. Acoust., Speech, Signal Process., vol. 33, no. 2, pp. 387–392, April 1985.
 [18] D. W. Rieken and D. R. Fuhrmann, “Generalizing MUSIC and MVDR for Multiple Noncoherent Arrays,” IEEE Trans. Signal Process., vol. 52, no. 9, pp. 2396–2406, Sep. 2004.
 [19] W. Suleiman, P. Parvazi, M. Pesavento, and A. M. Zoubir, “NonCoherent DirectionofArrival Estimation Using Partly Calibrated Arrays,” IEEE Trans. Signal Process., vol. 66, no. 21, pp. 5776–5788, Nov. 2018.
 [20] H. Kim, A. M. Haimovich, and Y. C. Eldar, “NonCoherent Direction of Arrival Estimation from MagnitudeOnly Measurements,” IEEE Signal Process. L, vol. 22, no. 7, pp. 925–929, Jul. 2015.
 [21] L. M. Kaplan, “Global node selection for localization in a distributed sensor network,” IEEE Trans. Aerosp. Electron. Syst., vol. 42, no. 1, pp. 113–135, Jan. 2006.
 [22] D. Li, K. D. Wong, Y. H. Hu, and A. M. Sayeed, “Detection, classification, and tracking of targets,” IEEE Signal Process. Mag., vol. 19, no. 2, pp. 17–29, Mar. 2002.
 [23] F. Elvander, I. Haasler, A. Jakobsson, and J. Karlsson, “Tracking and Sensor Fusion in Direction of Arrival Estimation Using Optimal Mass Transport,” in 26th European Signal Processing Conference, Rome, Italy, Sep. 37 2018.
 [24] T. T. Georgiou, J. Karlsson, and M. S. Takyar, “Metrics for power spectra: an axiomatic approach,” IEEE Trans. Signal, vol. 57, no. 3, pp. 859–867, Mar. 2009.
 [25] F. Elvander, A. Jakobsson, and J. Karlsson, “Interpolation and Extrapolation of Toeplitz Matrices via Optimal Mass Transport,” IEEE Trans. Signal. Process, vol. 66, no. 20, pp. 5285 – 5298, Oct. 2018.
 [26] F. Elvander, A. Jakobsson, and J. Karlsson, “Using Optimal Mass Transport for Tracking and Interpolation of Toeplitz Covariance Matrices,” in 43rd IEEE Int. Conf. on Acoustics, Speech, and Signal Processing, Calgary, Canada, April 1520 2018, pp. 4469–4473.
 [27] M. Cuturi, “Sinkhorn distances: Lightspeed computation of optimal transport,” Proc. Adv. Neural Inf. Process. Syst., pp. 2292–2300, 2013.
 [28] J. Karlsson and A. Ringh, “Generalized Sinkhorn iterations for regularizing inverse problems using optimal mass transport,” SIAM Journal on Imaging Sciences, vol. 10, no. 4, pp. 1935–1962, 2017.
 [29] T. T. Georgiou, “Solution of the general moment problem via a oneparameter imbedding,” IEEE Trans. Autom. Control, vol. 50, no. 6, pp. 811–826, 2005.
 [30] D. H. Johnson and D. E. Dudgeon, Array Signal Processing: Concepts and Techniques, Prentice Hall, Englewood Cliffs, N.J., 1993.
 [31] C. Villani, Optimal transport: old and new, Springer Science & Business Media, 2008.
 [32] JA. Luo, K. Yu, Z. Wang, and YH. Hu, “Passive source localization from array covariance matrices via joint sparse representations,” Neurocomputing, vol. 270, pp. 82–90, 2017.