Non-Coherent Sensor Fusion via Entropy Regularized Optimal Mass Transport

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 low-complexity iterations of the solution algorithm and thus makes the proposed method applicable also to higher-dimensional 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.


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 non-planar. Classical methods, originally used mainly for DoA estimation, such as beamforming techniques [5] and subspace methods like MU-SIC [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 well-calibrated arrays, i.e., known sensor positions, in order to produce reliable estimates. However, This work was supported in part by the Swedish Research Council, Carl Trygger's foundation, and the Royal Physiographic Society in Lund. 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, so-called noncoherent processing (see, e.g., [17][18][19][20]), i.e., we assume no information about cross-correlation 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, cross-correlations 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.

Near field localization
Consider an array consisting of p sensors located at y k ∈ R d , for k = 1, 2, . . . , p, and assume that point-like sources, located in a space S ⊂ R d , emit spherical waves impinging on the sensors. Then, the spatial distribution of the energy of the source signals may be representated by a non-negative function, or measure, Φ on S, referred to as the power spectrum. We will denote the space of such functions by M + (S). Letting λ denote the wavelength of the impinging waves, the covariance matrix of the sensor measurements is R = Γ(Φ), where the linear operator Γ : M(S) → C p×p is given by (see, e.g., [29]) where a(x) ∈ C p is the array manifold vector [30], i.e., with dm(x) denoting the Lebesgue measure. In other words, a positive semi-definite matrix R is a valid covariance for a sensor array if and only if R = Γ(Φ), for some Φ ∈ M + (S). From this, it may be noted that performing source localization based on sensor measurements corresponds to inferring a power spectrum Φ ∈ M + (S) from observations of the array covariance matrix R.
We will here consider a scenario wherein several sensor arrays measure the superposition of the impinging waves, i.e., when covariance matrices R j , for j = 1, . . . , J, are available, corresponding to J 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., Γ j (Φ) = R j , for j = 1, . . . , J, where Γ j is the linear operator corresponding to array j. However, for the case of miscalibration, the (assumed) operators Γ j contain errors, which may result in there being no single Φ ∈ M + (S) consistent with all covariance matrices. We will here adress this problem by utilizing OMT in order to induce robustness to the sensor fusion problem.

Optimal mass transport
Let Φ 0 and Φ 1 be two non-negative mass distributions on the space S. A transport plan M ∈ M + (S 2 ) for such a pair is a non-negative distribution on S 2 , with S 2 = S × S, where M (x 0 , x 1 ) describes the amount of mass transported from location x 0 ∈ S to location x 1 ∈ S. By defining a cost function c : S 2 → R describing the cost of moving mass, we may define the total cost associated with a transport plan as The OMT problem is then to find the feasible transport plan from Φ 0 to Φ 1 incurring minimal transportation cost (see, e.g., [31]), i.e., to solve where the constraints ensure that the transport plan is valid, i.e., matches the marginals Φ 0 and Φ 1 . 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 T (Φ 0 , Φ 1 ) will be used for interpolating between mass distributions Φ 0 and Φ 1 . It is worth noting that the OMT cost quantifies location errors based on the distances on the underlying space, in this case S ⊂ R d , thereby allowing for offsets due to, e.g., calibration errors. This is in contrast to standard metrics (e.g., L 1 , L 2 ) which only compare functions point by point without taking the geometry of the underlying space, i.e., relative location of the points, into account.

PROPOSED METHOD
Consider the localization scenario with the source signals impinging on a set of J arrays, with array manifold vectors a j (x), generating the covariance matrices R j , j = 1, . . . , J.
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 J 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 where Γ j denotes the linear operator (1) of array j. In order to find practical solutions to the barycenter problem in (5), the space S may be discretized as to be represented by n points x k ∈ S. Accordingly, the marginal spatial spectra are represented by the vectors Φ j ∈ R n + , and the cost function and transport plans are represented by matrices C = [c k, ] ∈ R n×n and M = [m k, ] ∈ R n×n + , respectively, where c k, = c(x k , x ), and m k, denotes the mass transported from x k to x . Although such a discretization yields finite-dimensional 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 where D(M ) = n k, =1 (m k, log(m k, ) − m k, + 1) is an entropy term, > 0 a small constant, and 1 ∈ R n denotes a vector of ones. In [27], it was shown that this problem may be solved in a computationally efficient way, using so-called Sinkhorn iterations.

Efficient sensor fusion
In order to arrive at a computationally efficient scheme for computing spectral barycenters, we propose to regularize the problem in (5)  subject to A j Φ j = ∆ j + r j for j = 1, . . . , J.
Here, r j = vec(R j ) ∈ C p 2 j , for j = 1, . . . , J, are the vectorizations of the covariances matrices, and A j ∈ C p 2 j ×n are matrix representations of the discretized operators Γ j . Also, the vectors ∆ j ∈ C p 2 j allow for noisy covariance matrix estimates, penalized by the standard 2 -norm, as controlled by the positive penalty parameter γ. The addition of the entropy term, D(M ), 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 1. For j = 1, . . . , J, let λ j be the solution to Algorithm 1 Sinkhorn-Newton method for the optimal mass transport barycenter problem with partial information. Require: Initial guess λ j , v j , and let u j = exp(A j λ j / ) for j = 1, . . . , J while Sinkhorn not converged do for j=1:J do while Newton not converged do

Update
v j = y./(K T u j (λ j )) for j = 1, . . . , J, Here, K = exp(−C/ ), and exp(·), ./, and denotes the elementwise exponential function, elementwise division, and elementwise product, respectively. The barycenter solution is then given by Φ = y, where y is the limit point of the iterations. Note that we here convert all complex quantities to an equivalent real-valued representation such that the vectorized covariances, r j , are mapped according to where Re(·) and Im(·) denote the real and imaginary parts, respectively, and where ∆ j and the discretized operators A j 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 noise-free case with complete observations of the underlying spectra, i.e., ∆ j = 0 and A j = I, as well as J = 2, one obtains the standard Sinkhorn iterations.  Fig. 1. Spectral estimate as given by the proposed estimator in (7). The alignment error is 6.7 degrees.
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.

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 [−0.5, 0.5] × [−0.5, 0.5]. 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 γ = 0.01 and = 0.005, and the cost function , for grid points x k , x ∈ R 2 . As comparison, we also obtain estimates using the noncoherent MUSIC and MVDR estimators, as described in [18], as well as the least-squares (LS) estimator from [32], and the non-coherent 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 devi-ations 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 non-perturbed 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.
A. PROOF OF PROPOSITION Proof 1 In the following, we denote the set {M 1 , . . . , M J } by M, and similarly for ∆, λ and µ. The Lagrangian of (7) with the dual variables λ and µ is Minimizing the Lagrangian with respect to M j gives the expression for the mass transport matrix solutions, with K = exp(−C/ ), u j = exp(A T j λ j / ) and v j = exp(µ j / ), for j = 1, . . . , J. Further, minimizing (9) with respect to the variables Φ and ∆, one gets a dual problem to where L(Φ * , M * , ∆ * , λ, µ) j = (λ T j r j − e (A T j λj ) T / e −C/ e µj / − 1 4γ λ T j λ j + n 2 ).
A block coordinate ascent method is to maximize (11) iteratively with respect to λ and µ. Maximization with respect to λ j requires r j − A j diag(e −C/ e µj / )e A T j λj / − 1 2γ λ j = 0, which gives the first step (8). To maximize (11) with respect to µ j , consider the corresponding Lagrangian with dual variable y, given by Maximization with respect to µ j requires −diag(e −C T / e A T j λj / )e µj / + y = 0 ⇔ v j = y./(K T u j ).
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