Scalable Gas Sensing, Mapping, and Path Planning via Decentralized Hilbert Maps.

This paper develops a decentralized approach to gas distribution mapping (GDM) and information-driven path planning for large-scale distributed sensing systems. Gas mapping is performed using a probabilistic representation known as a Hilbert map, which formulates the mapping problem as a multi-class classification task and uses kernel logistic regression to train a discriminative classifier online. A novel Hilbert map information fusion method is presented for rapidly merging the information from individual robot maps using limited data communication. A communication strategy that implements data fusion among many robots is also presented for the decentralized computation of GDMs. New entropy-based information-driven path-planning methods are developed and compared to existing approaches, such as particle swarm optimization (PSO) and random walks (RW). Numerical experiments conducted in simulated indoor and outdoor environments show that the information-driven approaches proposed in this paper far outperform other approaches, and avoid mutual collisions in real time.


Introduction
Fugitive emissions and the dispersion of pollutants in the atmosphere are significant concerns affecting public health as well as climate change. The accidental release of hazardous gases from both urban and industrial sources is responsible for a variety of respiratory illnesses and environmental concerns [1,2]. Many sensors are currently fabricated and deployed both indoors and outdoors for air quality data collection and communication. However, obtaining a spatial representation of a gas distribution is a challenging problem because existing mapping and fusion algorithms do not scale to networks of potentially hundreds of mobile and stationary sensors. The decentralized mapping and path-planning methods presented in this paper are applicable to very large-scale sensor systems, and as such can potentially be used to assess air quality, classify safe and hazardous areas based on the concentration of the harmful gases, and even localize fugitive emissions [3].
Due to the fundamental mechanisms of atmospheric gas dispersion, auxiliary tools are required to ensure early detection and to respond with appropriate counter measures by planning future measurements efficiently in both space and time. Methods for obtaining gas distribution maps (GDM) have recently been developed along two lines of research [4]. In the first line of research, a stationary sensor network is used to collect and fuse measurements to estimate the position of a source, typically requiring expensive and time-consuming calibration and data-recovery operations [5]. The use of

Problem Formulation
Consider the problem of optimally planning the trajectory of a distributed sensing system comprised of N cooperative robots engaged in GDM through a large region of interest (ROI), denoted by W ⊂ R 2 . In general, the gas concentration at a position x ∈ W, at time t, is modeled as a random variable C(x, t), thus the whole gas concentration is a random field. The gas distribution is considered constant over a time interval for simplicity and thus is modeled as a spatial function c(x) defined over W. The approach can be extended to time-varying concentrations by augmenting the dimensionality of the GDM.
Let U ⊂ R m denote the space of admissible actions or controls. The dynamics of each robot are governed by stochastic differential equations (SDEs), x n (t) = f[x n (t), u n (t), t] + Gw(t), x n (T 0 ) = x n 0 and n = 1, · · · , N where x n (t) ∈ W denotes the nth robot's position and the velocity at time t, u n (t) ∈ U denotes the nth robot action or control, and x n 0 denotes the robot initial conditions at initial time T 0 . The robot dynamics in (1) are characterized by additive Gaussian white noise, denoted by w(t) ∈ R 2 , and G ∈ R 2×2 is a known constant matrix. In this paper, we assume that the position of the nth robot is obtained by an on-board GPS and is denoted byx n (t), where (·) denotes the estimated variable's value. All N robots are equipped with identical metal oxide sensors in order to cooperatively map a time-constant gas distribution, where the quantity and position of gas plumes are unknown a priori. Because the reaction surface of a single gas sensor is very small (≈1 cm 2 ), a single measurement from a gas sensor can only provide information about a very small area. Therefore, to increase the resolution of the GDM, a small metal oxide sensor array is used instead of a single oxide sensor for each robot. The area that is covered by this small metal oxide sensor array is called the field of view (FOV) of the robot, denoted by S(x n ) ∈ W. Because the structure of the metal oxide sensor array is fixed, the positions of each metal oxide sensor can be approximated by the robot position measurements, denoted byx n,m , m = 1 : M, where M is the number of oxide sensors on-board the nth robot. For example, for a 5 × 5 metal oxide sensor array, the number of sensors is M = 25.
The gas concentration measurement obtained by the mth sensor can be modeled as the sum of the actual concentration and measurement noise, c n,m (x n ) = c(x n,m ) +v(x n,m ) (2) wherev(x n,m ) is a realization of the random measurement noise, V(x n,m ). Furthermore, the concentration measurement,ĉ n,m (x n ), can be treated as a sample of the random variable C(x n,m ), defined as Then, the concentration measurementĉ n,m (x n ) can be considered as an observation of the random field C(·) at position x n,m .
Because the actual gas concentration, c(x), and the distribution of the measurement noise, V(x), are unknown, our objective is to approximate the random field C(·), and then use this approximated random field to estimate the gas concentration map. The cost function of the nth robot over a fixed time interval [t 0 , t f ] can then be expressed as an integral of the future measurement values, conditioned on past measurements, and the robot control usage, i.e., where H[C n (x|M n (t))] is the entropy of the random variable C n (x|M n (t)) approximated by the nth robot given all past measurements, M n . R is a positive definite matrix that weighs the importance of the elements of the control input u n (t), and the superscript "T" denotes the matrix transpose.
The optimal planning problem is to find the optimal time history of the robot state x n (t) and control u n (t), for all n = 1, · · · , N, so that the cost function J n in (4) is minimized over [t 0 , t f ], and subject to (1). Let the time interval [t 0 , t f ] be discretized into T f equal time steps ∆t = (t f − t 0 )/T and let t k = t 0 + k∆t denote the discrete-time index with k = 1, · · · , T f . Then the objective function, J n , can be rewritten as, where the subscript k indicates the kth time step and the superscript "(T f )" indicates the time until the T f th step. The term, M M n,k , indicates the measurement history until the T f th step and M n,k denotes all the measurements obtained by the nth robot up to the kth time step.

Representation of GDM
To approximate the probability distribution of the continuous random variable C n (x|M (T) n ), there are several methods, including kernel density estimation (KDE) [24,25] and Gaussian process regression (GPR) [26]. In the KDE method, the learned probability density function (PDF) is assumed to be a weighted summation of many parameterized kernel functions, where the weight coefficients and the kernel parameters are critical and learned from the raw measurement data set. In the GPR method, the approximated PDF is assumed to be a Gaussian distribution, where a large matrix, the Gram matrix, is very critical and learned from the raw measurement data. For both of methods, the data fusion of two different measurement data sets and update of the approximated PDF are computationally demanding, because all computations must be implemented from a massive amount of raw data. It is difficult to obtain the updated parameters and coefficients from the previous parameters and coefficients directly.
The method in this paper overcomes this hurdle through an efficient fusion approach developed by approximating the continuous probability distribution by a discrete probability distribution. Denote the range of the concentration in the whole ROI by R = [L c , H c ], and then divide the range into L concentration intervals denoted by R 1 = [L 1 , H 1 ), · · · , R l = [L l , H l ), · · · , R L = [L L , H L ], where H l = L l+1 for l = 1, · · · , L − 1, and L 1 = L c and H L = H c . The cutoff coefficients, {(L l , H l )} L l=1 , specify the concentration intervals of interest. Instead of approximating the PDF, this paper models the probabilities that the gas concentration at every position x ∈ W belongs to the concentration interval {p l (x)} L l=1 , or, where P(·) denotes the probability operator. Let f c (x) denote the PDF of the concentration C(x) and ∆L l = H l − L l denote the length of the lth concentration interval. The relationship between p l (x) and f c (x) can be expressed as, where Therefore, the discrete distribution {p l (x)} L l=1 can be used to describe the probability distribution of the continuous random variable C(x) as L → ∞. More importantly, it can be used for decision making, where the concentration intervals correspond to different levels of hazardous gas concentration. The Hilbert map method, developed by Fabio Ramos and Lionel Ott in [27] to model obstacle occupancy, is modified here to approximate the discrete distribution {p l (x)} L l=1 over the entire ROI.

Mapping with KLR
A Hilbert map is a continuous probability map developed by formulating the mapping problem as a binary classification task. Let x ∈ W be any point in W and Y ∈ {0, 1} be defined as a categorical random variable such that where the realization is denoted by y. The Hilbert map describes the probabilities P(Y = 1|x) = p(x) and P(Y = 0|x) = 1 − p(x) at the position x.
Consider the training measurement data set, M = {(x m , y m )} M m=1 , where x m ∈ W indicates the measurement position, the Boolean variable y m ∈ {1, 0} indicates whether the concentration measurement,ĉ(x m ), belongs to the range, R (where 1 = Yes and 0 = No), and M is the number of measurements. The probability P(Y|x) is defined as, where, (11) and f ∈ H is a Hilbert function defined on W. H is a reproducing kernel Hilbert space (RKHS) associated with a kernel k(·, ·). The kernel mapping is denoted by k(x, ·) = ϕ(x), ϕ(x) ∈ H, and is an injective map from W to H [28][29][30][31][32]. According to the kernel trick [28][29][30][31][32], the evaluation of the Hilbert function can be expressed in the form of inner product, where ·, · H indicates the inner product in the RKHS. To learn the Hilbert map, the loss function J H is defined as, where m = − ln [P(Y = y m |x m )] is the negative log-likelihood (NLL) of the data (x m , y m ). Here, the term λ is the regularization term, which is a small user-defined positive scalar, λ 1. Then the gradient of the loss function with respective to f is expressed as, T , and y = [y 1 , · · · , y M ] T . In addition, the Hessian operator is expressed as, where I is an identity operator (or matrix) defined in the domain of H × H, so that for any function h ∈ H, and Ih = h. Here, Λ is an M × M diagonal matrix defined as, and 1 M = [1, · · · 1] T is an M × 1 vector with all elements equal to one. Using the Newton-Raphson method, the Hilbert function, is updated iteratively, where where It can be seen from (18) that the evaluation of f i+1 (x) only depends on the evaluations of f i (x) and where f 0 (x) is the initial function evaluation at x, prior to learning the function from M. The following matrix only depends on the measurement data, M, and can be updated iteratively. Furthermore, consider Q collocation points, X c = {x c q ∈ W } Q q=1 , in the ROI, characterized by the same spatial interval, labeled by the superscript "c". Let Φ c = ϕ(x c 1 ), · · · , ϕ(x c Q ) denote the feature matrix of the collocation points. The evaluations of function f (x c q ) at all collocation points can be updated by, where According to (19), the evaluations, where f 0 comprises the initial function evaluations at the collocation points prior to learning the function from the measurement data M.
In summary, instead of learning the function f or its coefficients as in traditional KLR or GPR [26], we update the evaluations of f (x) at the collocation points, X c . Given the Hilbert function f , the evaluations at the collocation points provide the Hilbert map, f, defined as

Temporal Update of Hilbert Map
The previous subsection develops a method for learning the Hilbert map from the measurement data set M, obtained by the sensors during one time interval. This subsection considers a Hilbert map updated according to the next data set, where λ is the regularization term as in (13), γ(T − k) is the "forgetting factor", and k,m = − ln [P(Y = y k,m |x k,m )] is the NLL of the data (x k,m , y k,m ). Similarly to (14), the gradient is expressed as, on the diagonal of a zero matrix. In addition, as with (15), the Hessian operator can be expressed as, where,Λ and Then, the update rule for learning the Hilbert function, f , at the ith iteration is given by, According to the above equation, the function, f , is expressed by using the set of all the measurement positions, { T k=1 X k }, and the corresponding coefficients.
To reduce the computational load, the forgetting factor, γ(T − k), is modeled by, such that each robot stores only the measurements obtained during the past τ time steps. By setting τ = 1, the update rule (28) for learning the function f at the ith iteration is expressed as, where

Approximation of GDM
The previous subsections present a method for learning the Hilbert function, f , from the available measurement data, comprising the measurement position data set X k = {x k,m } M k m=1 and corresponding classification estimates Y k = {y k,m } M k m=1 for k = 1, · · · , T that indicate whether the concentration measurement at the measurement position belongs to the range, R, defined in (9). If L classification estimates are obtained from the same concentration measurements, indicating that the concentration measurement belongs to the L concentration ranges, such as R 1 , · · · , R L , respectively, using the learning method described in Section 3.2, one can approximate L Hilbert functions, f 1 , · · · , f L . Furthermore, the probabilities in (6) can be evaluated at all the collocation points, X c = {x q } Q q=1 , such that p l (x c q ) = 1 − 1 e f l (x c q ) , q = 1, · · · , Q and l = 1, · · · , L Now, consider Hilbert functions, { f l } L l=1 , approximated locally by each robot. Although all the concentration intervals are mutually exclusive and complete, i.e., it is not guaranteed that the sum of all the learned probabilities, {p l (x c q )} L l=1 , is equal to one, or ∑ L l=1 p l (x c q ) = 1. Therefore, the learned probabilities must be normalized to obtain the discrete distribution, π(x c q ) = [π 1 (x c q ), · · · , π L (x c q )]. Each component of the discrete distribution is calculated by, where p 0 (x c q ) is a normalization term. Letc = [c 1 , · · · ,c L ] denote the medians of these intervals, wherec l = (H l − L l )/2. Then, the expectation of the random variable C(x c q ) can be expressed as, where E(·) is the expectation operator.

Information Fusion
In this section, the problem of information fusion between neighboring robots is considered, where each robot builds its own Hilbert function locally using a decentralized approach. Assume that all the robots use the same collocation points, X c . To limit the amount of communicated data, the evaluations of the Hilbert functions at the collocation points are shared among the robots, instead of the coefficients and parameters of the Hilbert functions.

Hilbert Map Fusion
Consider two robots that have each learned their respective Hilbert functions, f 1 and f 2 , from two different sets of measurement data, M 1 and M 2 , respectively. Then, the fused Hilbert function, f F , can be obtained based on the following theorem.
Theorem 1. Let f 1 (x) and f 2 (x) be two Hilbert functions defined on a workspace W, and approximated by two robots based on their own measurement data sets, M 1 and M 2 , respectively. These two Hilbert functions can be applied to calculate the conditional probability that the concentration is in the range R given the corresponding measurement data sets, as follows: Then, the fused conditional probability p F (x|M 1 , M 2 ) can be expressed as where f F (x) is the fused Hilbert function. In addition, the fused Hilbert function, f F (x), can be calculated from the Hilbert functions, f 1 (x) and f 2 (x), as follows, is the ratio between the prior probabilities, P(C(x) ∈ R) = p(x) and P(C(x) / ∈ R) = 1 − p(x), at x ∈ W.
The proof of Theorem 1 is provided in the Appendix A.1. According to Theorem 1, the following corollary can be obtained.

Corollary 1.
Assume that the prior probability is even, such as p(x) = 1/2. Then, the ratio between the prior probabilities can be calculated by, and the fusion function can be rewritten as Because the prior concentration distribution is unknown, the Hilbert functions are fused according to (41) in Corollary 1, unless otherwise stated. Furthermore, assume that all the robots have the same collocation points, X c . The information fusion can be implemented by fusing the Hilbert maps among neighboring robots, where f H,1 and f H,2 , and f H,F are the Hilbert maps associated with the Hilbert functions f 1 (x), f 2 (x), and f F (x), respectively.

Communication Strategy
Any pair of robots in the network can share and update their Hilbert maps efficiently according to (42). To implement data fusion for a very large number of robots, however, a communication strategy is also required to determine how to share data between robots characterized by active communication links. In this paper, gas measurement data are communicated at every time step T c . The communication protocol requires four steps at every time k, as follows. In the first, the N robots are partitioned into smaller communication networks according to their positions and communication range r c , and G ı denotes the index set of the robot in the ıth communication network. Then, for any robot, a n , with n ∈ G ı , there exists another robot, a n , such that where d n,n is the distance between robots a n and a n , and · indicates the Euclidean norm. As a second step, one robot in every communication network is selected as a temporary data center (TDC) denoted by a n * ı , n * ı ∈ G ı . The other robots in G ı send the Hilbert map change ∆f n,k to robot a n * ı . The Hilbert map change, ∆f n,k , is defined as the change between the nth robot's Hilbert map at the current time step k and its Hilbert map at the previous communication time step, k − T c , such that where f n,k and f n,k−T c denote the Hilbert maps of the nth robot at times kth and (k − T c )th, respectively. In the third step, the sum of the Hilbert map changes obtained from all robots in G ı , defined as, is communicated back to the other robots, except the TDC robot, a n * ı . Finally, in the fourth step, all the robots update their own Hilbert maps by adding the received total Hilbert map changes, ∆f G ı ,k , to the current Hilbert maps, such that f n,k + = f n,k + ∆f G ı ,k , n ∈ G ı where f n,k + represents the Hilbert map after the data fusion process.

Path-Planning Algorithms
In the previous sections, the Hilbert maps f l,n,k , l = 1, .., L, n = 1, · · · , N and k = 1, · · · , T f , are approximated and updated by the robots. The corresponding binary probabilities p l,n,k can be calculated efficiently as follows p l,n,k = p l,n,k (x c 1 ), · · · , p l,n,k (x c Q ) and the entropy at each collocation point x c q ∈ X c , q = 1, · · · , Q, is obtained by Then, an entropy map, h n,k is obtained from the vector, where h n,k (x c q ) = ∑ L l=1 h l,n,k (x c q ) denotes the sum of the entropy at the collocation point x c q . According to the cost function in (5), the objective of the nth robot is to minimize the uncertainty of the concentration distribution, which can be implemented by minimizing the entropy at all the collocation points, such thatJ Therefore,J n (T f ) can be treated as an approximation of the term W H[C n (x|M (5). In the rest of paper,J n (T f ) is applied as the approximation of the final cost function in (5) unless otherwise stated. In other words, the nth robot should visit the area around the collocation point x c q , which has the higher value of h n,k (x c q ) at the kth time step. Based on this idea, two entropy-based path-planning algorithms are proposed in the following section to control the robots such that the value of the concentration measurements obtained by the robots in the ROI can be optimized over time.

Entropy-Based Artificial Potential Field
An information-driven approach is developed by planning the path of the robots such that they move towards collocation points with higher entropy. The collocation points are treated as goal positions characterized by attractive artificial potential fields defined as, where the superscript "att" indicates the attractive field. The corresponding attractive gradient is expressed as Similarly to classic artificial potential field methods, the repulsive potential functions U rep generated from the other robots are also considered, such that where r rep is the distance threshold to create a repulsion effect on the robot. The repulsive gradient is expressed as Using (52) and (54), the total potential gradient for the nth robot is expressed as, where ξ and η are user-defined coefficients. Based on the total gradient, g tol (x n ), the nth robot can be controlled to visit the collocation points with higher uncertainty and avoid collisions with other robots. The algorithm is developed based on the entropy attraction, which is named as EAPF algorithm.

Entropy-Based Particle Swarm Optimization
The particle swarm optimization (PSO) algorithm proposed by Clerc and Kennedy in [33] and its variants use a "constriction coefficient" to prevent the "explosion behavior" of the particles, and have been successfully applied to GDM and gas source localization problems [12,13]. In the original PSO algorithm and its variants, the concentration measurements are used to update the robot controls. Considering the different objective function in (5), an entropy-based PSO (EPSO) is proposed.

Simulations and Results
The performance of the decentralized GDM methods presented in this paper is demonstrated on a gas sensing application, where information about gas concentration obtained by a large network of robots is fused and used for information-driven path planning in a decentralized approach. Two indoor and outdoor environments are simulated and used to test the proposed methods. The new entropy-based EAPF and EPSO path-planning algorithms are compared to the existing algorithms known as random walk (RW) and classical particle swarm optimization (CPSO) [12,13,33].

Indoor GDM Sensing
The decentralized sensing system consists of a network of N = 100 robots characterized by single integrator dynamics,ẋ where x i = [x, y] T is the robot state vector, x, y are the inertial coordinates, and u i = [u 1 , u 2 ] T is the robot control vector comprised of the xand y-velocity components. The robot state/position is assumed to be observable by a built-in GPS with zero-mean measure noise w, where w is Gaussian white noise N (0, Σ) with Σ = 0.05 × I 2 . The above  Figure 1. Each robot is equipped with a small metal oxide sensor array comprised by M = 5 × 5 = 25 gas sensors, where the spatial intervals between two sensors are all 10 cm. The FOV of each sensor array covers an area of 40 × 40 cm 2 in the ROI, which is very small relative to the whole ROI. Assume that the measurement noise of the gas concentration V(x) is also characterized by Gaussian white noise, N (0, Σ c (x)) with the covariance matrix Σ c (x) = 0.5 × I 2 everywhere in W. All of robots communicate with each other at the same time with a communication period T c = 10 min. The communication range is r c = 30 m. In addition, Q = 100 × 100 virtual collocation points are evenly deployed in the ROI to generate the Hilbert maps, where the Gaussian kernel is used for Hilbert function learning with a kernel size of σ = 2 m, and the parameter τ is set to 3. The EPSO and the CPSO neighbor range coefficient, r neig = 5 m, is applied to determine the best neighboring collocation points. The maximum velocity of each robot is 2 m/min for all simulations. The four path-planning algorithms, referred to as EAPF, EPSO, CPSO, and RW, are tested for mapping the gas distribution. The approximated cost functionJ n (k) defined in (50) is calculated by each robot at every time step as shown in Figure 2, where the solid line represents the mean of the approximated cost values over all of the robots at the kth time step, and the dashed line indicates one standard deviation above and below the mean. These cost histories reflect how effective the robots are at mapping the gas distribution. A lower value ofJ n (k) indicates that more information about the gas distribution is obtained by the robot at the kth step. It can be observed that the EAPF and EPSO algorithms achieve significantly better performance than the others. They can map the gas distribution more rapidly and completely, while the CPSO and RW algorithms perform poorly because they do not used prior measurement data for planning.
The means and standard deviations of the approximated cost function values of the all robots at the final time k = T f = 500 for all the algorithms are tabulated in Table 1. It can be seen that the approximated cost function obtained by the EPSO algorithm outperforms the other algorithms in the indoor environment. Let n * denote the index of the robot who gets the lowest value ofJ n (T f ) at the final step. Then, the value ofJ n * (k) represents the best performance among all the robots. The estimates of the gas concentration in the ROI are calculated based on these learned Hilbert maps according to (35), wherec = [15,50,85] is calculated from the cutoff coefficients. For each path-planning algorithm, three snapshots of the estimated GDMs generated by the n * th robot in each simulation are presented in Figure 3, where the snapshots are taken at the k = 20th, k = 100th, k = 500th time steps, respectively. It can be seen that the robots controlled by the EAPF and EPSO algorithms obtain clean GDMs at the final time step, while the robots controlled by the other two existing algorithms cannot complete the mapping task in the given time period.  The normalized mean square errors (NMSE) between the estimated gas distribution and the actual gas distribution are calculated for each robot. The means and the corresponding standard deviations of the NMSE over the all robots for the different planning algorithms are reported in Table 2, which obviously shows that the EPSO algorithm outperforms the other algorithms to estimate the GDMs.

Outdoor GDM
To verify the robustness and versatility of the proposed approaches, a GDM shown in Figure 4, originally presented in [34], is used to represent the gas distribution in an outdoor environment. The gas concentrations are normalized to the range R = [0, 100] for comparison like the indoor simulations. In this case, the intervals are chosen as R 1 = [0, 60), R 2 = [60, 80), and R 3 = [80, 100], to represent the plume shapes. All other parameters, including the initial robot positions, are the same as those used in the previous subsection.
The approximated cost functionJ n (k) obtained by all four algorithms is plotted in Figure 5. The approximated cost function at the final step are also tabulated in Table 3. Similarly to the indoor simulations, the gas concentration is estimated according to (35), wherec = [30,70,90] is calculated based on the cutoff coefficients. The evolution of the estimated GDM obtained by different algorithms is presented in Figure 6. Furthermore, the statistical results of the NMSE between the estimated GDM and the actual GDM are reported in Table 4.
x (m) y (m) Normalized gas concentration     As expected, the results presented in Figures 5 and 6 and Tables 3 and 4 all show that the proposed EAPF and EPSO algorithms work well in the outdoor GDM problem, while the CPSO and RW algorithms cannot map the gas concentration in the entire workspace in the given time period. In addition, the EPSO algorithm significantly outperforms the other algorithms in all simulations.

Conclusions
This paper presents a decentralized framework for GDM and information-driven path planning in distributed sensing systems. GDM is performed using a probabilistic representation known as a Hilbert map and a novel Hilbert map fusion method is presented that quickly and efficiently combines information from many neighboring robots. In addition, two entropy-based path-planning algorithms, namely the EAPF and EPSO algorithms, are proposed to efficiently control all the robots to obtain the gas concentration measurements in the ROI. The proposed approaches are demonstrated on a system with hundreds of robots that must map a gas distribution collaboratively over a large ROI using on-board iron-oxide arrays and no prior information. The results show that through fusion and decentralized processing, the entropy of the gas map decreases over time, the robot paths remain safe (avoiding mutual collisions), and the entropy-based methods far outperform both traditional and random approaches.

Conflicts of Interest:
The authors declare no conflict of interest and the funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Given two Hilbert functions f 1 (x) and f 2 (x), the conditional probabilities that the concentration at the position x is belong the range R are expressed as p 1 (x|M 1 ) = P(C(x) ∈ R|M 1 ) (A1) p 2 (x|M 2 ) = P(C(x) ∈ R|M 2 ) (A2) For convenience, the event C(x) ∈ R and its complement even C(x) / ∈ R are denoted by e and e C , respectively. Then, the fused conditional probability, p F (x|M 1 , M 2 ), can be expressed as, = P(e)P(M 1 |e)P(M 2 |e) P(e)P(M 1 |e)P(M 2 |e)+P(e C )P(M 1 |e C )P(M 2 |e C ) = P(e)P(M 1 |e)P(e)P(M 2 |e) P(e) 2 P(M 1 |e)P(M 2 |e)+εP(e C ) 2 P(M 1 |e C )P(M 2 |e C ) = P(e,M 1 )P(e,M 2 ) P(M 1 )P(e|M 1 )P(M 2 )P(e|M 2 )+εP(M 1 )P(e C |M 1 )P(M 2 )P(e C |M 2 ) = P(e|M 1 )P(e|M 2 ) P(e|M 1 )P(e|M 2 )+εP(e C |M 1 )P(e C |M 2 ) = p 1 (x|M 1 )p 2 (x|M 2 ) p 1 (x|M 1 )p 2 (x|M 2 )+ε[1−p 1 (x|M 1