Structurally aware 3D gas distribution mapping using belief propagation: a real-time algorithm for robotic deployment

This paper proposes a new 3D gas distribution mapping technique based on the local message passing of Gaussian belief propagation that is capable of resolving in real time, concentration estimates in 3D space whilst accounting for the obstacle information within the scenario, the first of its kind in the literature. The gas mapping problem is formulated as a 3D factor graph of Gaussian potentials, the connections of which are conditioned on local occupancy values. The Gaussian belief propagation framework is introduced as the solver and a new hybrid message scheduler is introduced to increase the rate of convergence. The factor graph problem is then redesigned as a dynamically expanding inference task, coupling the information of consecutive gas measurements with local spatial structure obtained by the robot. The proposed algorithm is compared to the state of the art methods in 2D and 3D simulations and is found to resolve distribution maps orders of magnitude quicker than typical direct solvers. The proposed framework is then deployed for the first time onboard a ground robot in a 3D mapping and exploration task. The system is shown to be able to resolve multiple sensor inputs and output high resolution 3D gas distribution maps in a GPS denied cluttered scenario in real time. This online inference of complicated plume structures provides a new layer of contextual information over its 2D counterparts and enables autonomous systems to take advantage of real time estimates to inform potential next best sampling locations.


I. INTRODUCTION
A UTONOMOUS robots have seen a dramatic increase in research in recent years due to their ability to perceive environments without human intervention.In this regard, the use of mobile robotics within chemical, biological, radiological and nuclear (CBRN) applications is of great modern interest.Recent disaster events, such as the Fukushima disaster in 2011, and the ever present threat of chemical warfare have dictated the need for quickly deployable systems that can identify the characteristics of a gas dispersion event whilst operating in potentially unknown environments.In such scenarios, key media inferred by mobile systems can inform first responders whilst simultaneously providing safe data collection.
Gas distribution mapping answers the proposed question of How is the gas distributed across an area?Whilst in certain scenarios it may only be of interest to determine if a gas is present or where the source of the release is, the gas distribution across the whole environment can provide richer information.Contextually, distribution maps be can used to identify areas of high concentration (as well as low concentration) in order to dictate exclusions zones, the boundary of the plume and likely locations of a leak.This information is a key first step in hazard response.Furthermore, inferring distribution maps in real time as data are collected can help to determine the next best sampling location for a sensor in both the autonomous and teleoperated cases.
A combination of small scale fluctuations in the dispersal pattern, a limited sensing area and transient measurements associated with chemical sensors makes gas mapping an exceptionally difficult task compared to the similar CBRN related field of radiation mapping, which can rely on a smooth inversesquared relationship to make predictions [1].Equally, more complicated model based methods that utilise computational fluid dynamics (CFD) (e.g.[2]) are capable of producing high fidelity distributions that can account for complex geometry.However, these methods rely on reliable prior meteorological information and can take in the order of days to converge to an accurate solution.This makes them unsuitable for both first response and mobile sensor data collection which both require solutions in real time.Based on these insights, there has been an increasing trend towards investigating data driven methods for mapping a gas distribution.These methods make no strong assumptions on the gas transport model and instead rely on online fitting of a function to dense and diverse data.
Kernel DM, first proposed by Lilienthal in 2003 [3], is one of the first attempts to tackle the problem of concentration mapping for mobile robotics.This is achieved by arXiv:2212.07552v1[cs.RO] 14 Dec 2022 discretising a spatial map into a grid array with each cell in the grid storing the mean concentration predicted by the Kernel algorithm.Each measurement then applies a Gaussian shaped contribution to the surrounding cells based on the distance away from the measurement.By formulating the problem in this manner, concentrations in the environment are only dependent on the measurement themselves and not their surrounding cells.The work has been extended in [4], where Kernel DM+V provides a distribution of the confidence and the variance, where variance relates to the local variability of concentrations (source term predictor) and confidence relates to the uncertainty of the mean distribution.In [5], Kernel DM+V/W is developed.In this work, the Gaussian kernel, is transformed into a bivariate Gaussian along the axis of a measured wind direction, thus enabling the inclusion of anemometry sensor measurements.A further addition to the kernel method is shown in [6], where TD kernel DM+V/W allows for the time decay of measurements as they become historically uncertain.These works outline the core of the kernel DM algorithm, however further works have investigated multi-compound discrimination [11], a 3D kernel [12] and integrating obstacle effects [13].The obstacle inclusion works by scaling the contribution from a measurement to a cell by the Euclidean distance transform over the map from the measurement to the cell.This must be calculated for each measurement and also each cell within its radius and therefore is computationally inefficient for high resolution maps or 3D mapping tasks.Whilst Kernel DM is a lightweight local GDM algorithm, the method does not implicitly include the effects that obstacles have on the plume in a statistical manner due to the fundamental lack of dependency between predictions and the spatial map.
The GP technique as a general regression model, has also been applied to GDM in early work by Stachniss et.al [9] and more recently by Hutchinson et.al [10].In the earlier example, a Gaussian Process mixture model is employed wherein a mixture of 2 Gaussian processes are simultaneously learnt with one for the smooth background concentration function and one to model the plume concentration.The Matern 3/2 covariance function is used for the shaping kernel.The mixture model is compared against the Kernel DM method and single GP technique in an experimental trial and the performance of mixture model is shown to be the best at predicting the gas distribution.GP models are capable of predicting continuous distributions over a domain, contrary to discrete methods, and therefore are not susceptible to discretisation errors.However, GPs scale poorly with an increasing number of measurements (which is necessary for large scale mapping tasks) and are not currently capable of accounting for obstacle integration with the predicted map.
Most recently, the GMRF algorithm for GDM has been proposed by Monroy et.al [7] that discretises the environment into cells and then forms a graphical model that describes the relationship between adjacent cells and measurements.In fact, it is shown in that a GMRF is equivalent to a GP in [14].Contrary to the GP, due to modelling the dependency between a cell and its neighbours, the GMRF model probabilistically accounts for the effects of obstacles on the concentration at a given cell whilst also accounting for the increasing uncertainty of measurements over time (similarly to TD Kernel DM).The underlying graphical model used is known as a factor graph and is a common representation of inference problems in mobile robotics [15].Due to the Gaussian assumptions placed upon the factor graph, the problem reduces to a GMRF that can be solved by the least squares method.In [8], the algorithm is extended to also model the wind field via measured anemometer readings and resolve the effects that a wind field has on a cells gas concentration (bringing parity with TD Kernel DM+V/W).The major downside of G-GMRF stated in [7], is that the computational complexity is O(n 1.5 ) and the entire map must be resolved per iteration.This limits the size of map that a mobile robot can calculate and also requires that the map is fully resolved each time an update to the graph is added.
To overcome this issue, Rhodes et.al [16] advocate the use of the Gaussian belief propagation (GaBP) solver on the GDM factor graph problem.Belief propagation is a message passing technique that can be applied to loopy graph structures [17] (such as those factor graphs in GDM) and has seen increasing prominence in mobile robotics in recent years, especially in the SLAM domain [18]- [20].The G-GaBP algorithm proposed makes use of the wildfire message schedule [18] to iteratively solve the factor graph via local message passing propagating around sensor measurements.This shift away from global solvers (e.g., the direct method used in [7]) towards local inference allows for orders of magnitude quicker solutions to sequential mapping, and can also be used to acquire anytime estimates of the distribution maps (as opposed to waiting for a batch solution).Whilst the initial work on G-GaBP proposes the basic application of belief propagation to the GDM domain, simulation work is only carried out on a 2D scenario and inference is still performed on a globally defined map.
Therefore, we extend our work by addressing the problem of performing inference on a 3D gas distribution map, applying several efficiency upgrades, and also showing the real-world applicability of the algorithm to an online mapping robotic platform.We now list the specific contributions of the work displayed in this paper.

B. Contributions
This paper consists of four main contributions to the state of the art.Firstly, we formulate the GDM problem as a 3D network of connected states and measurements from which we can perform probabilistic inference, increasing functionality from the 2D state of the art.
We then propose two major improvements to the GaBP solver proposed in [16] to increase its efficiency when applied to mapping with multiple point sampling sensors.The first improvement concerns the schedule for sending messages in the belief propagation algorithm.A hybrid solver is proposed that utilises the properties of both the wildfire and residual belief propagation algorithms to increase the global convergence rate.The second improvement reformulates the inference problem from a static a-priori graph into a dynamic information theoretic graph structure that grows relative to the information inferred.
Our final contribution is the application of the proposed solver into a robotic framework that is capable of performing real-time inference in completely a-prior unknown environments.By doing so, we show the first mobile robot capable of performing 3D gas distribution mapping in conjunction with a dynamic occupancy map (to inform factor graph construction), whilst also performing real time simultaneous localisation and mapping (SLAM).All the above functions are shown being performed entirely onboard the robotic platform with real-time inference updates available to a potential operator or informative path planner.We also make the deployed algorithm available as on open source ROS package for others to implement on their systems.
The remainder of the paper is as follows.Section II formulates the GDM problem as a 3D factor graph and shows how it can be simplified to a GMRF under Gaussian assumptions.We then introduce the governing equations behind the GaBP framework in Section III and describe how the GMRF problem can be solved using local message passing.Next, the proposed 3D G-GaBP algorithm is derived in Section IV detailing the aforementioned developments to the propagation framework.The proposed algorithm is then tested in several simulations in Section V for quantitative analysis.Finally, the experimental evaluation is described in Section VI, the proposed algorithm is deployed onboard a robotic platform in a real-world 3D mapping task.

II. PROBLEM FORMULATION
Let X = {x 1 , . . ., x n } be a finite set of random variables that we wish to infer with each variable x i representing a distribution of the concentration at the i-th cell or voxel in the domain D ⊂ R 3 .We assume that the sensing robot takes a set of measurements at different locations within the domain.Denoting z k,i the noisy sensor measurement taken at time k associated with x i , the set of measurements corresponding to x i collected at different timestamps is denoted as Z i .With a slight abuse of notation, we use the subscript k as the index of the measurement within the set Z i .The total measurement set collected so far is expressed as Z = {Z i | Z i = ∅, i = 1, . . ., n}.We aim to resolve the maximum a posteriori (MAP) estimate of the variables x = x 1 • • • x n T from the joint probability density function p(x|Z) to give the most likely cell concentrations across the gas mapping domain.

A. Factor graphs
To model the dependencies between these random variables and the measurements recorded, the GDM problem can be formulated as an inference problem over a factor graph G = (X , F, ε), which is an undirected bi-partite graph containing variable nodes X , factor nodes F = {f 1 , f 2 , . . ., f S } and their edges ε encoding the connectivity according to the variables involved in each factor.Each factor f s (•) describes the dependencies between a local subset of variable nodes X s ⊂ X , which are directly connected to the factor.Considering all nodes in the graph, the joint probability distribution is expressed as the product of all factors in the graph [15]: Factors in the general sense can take any form, however, it is typical in mobile robotics to represent factors as Gaussian distributions, since it is a common assumption that sensor measurements deviate from the ground truth value in a Gaussian manner and that the gas concentration values also affect their surroundings spatially in a Gaussian manner, as in other GDM algorithms (e.g.Kernel DM and GP).When applying the Gaussian assumption to all factors in the factor graph, it can now be represented as a GMRF.From this model, the problem of retrieving the MAP solution of vector x is the same as finding the marginal means of the associated Gaussian joint probability distribution being p(x|Z) = N (Λ −1 g, Λ −1 ), i.e., where Λ is the Information matrix of the GMRF model and g is the observation vector [15].The construction of such a model for GDM is derived in section II-B.By factorising the pairwise cliques defined by the factors in the graph between variable x i and x j , Eq. ( 1) can be rewritten in GMRF form as a function of edge potentials, ψ ij exp(−x i Λ ij x j ), and self potentials,

B. Factor graph construction for 3D-GDM
In our previous work [16], we make use of the typical 2D graph construction of Monroy et al. [7].We now extend the factor graph formulation to account for 3D gas distribution maps.For the basic GDM problem, we take three types of Gaussian factors that represent the correlations between concentration variables and some observed quantity, including the relationship between neighbouring discrete cells in R 3 space.These are observation factors f o , regularisation factors f r and default factors f d .
The observation factor defines the potential between a noisy sensor measurement z k,i and the modelled concentration value x i .Since the measurements are observed variables and not part of the inference itself, these potentials contribute to a nodes self potential φ i and are often referred to as unary factors.The sensor measurement is modelled as z k,i = x i + ω k + ζ k with two additive Gaussian noises, of which ω k ∼ N (0, σ 2 s ) arises from the sensor noise and ζ k ∼ N (0, σ 2 ζ ∆t) is from the time dependent corruption that models the increase in uncertainty over time.σ 2 ζ should be set relative to the volatility of the underlying gas distribution (i.e. higher for more volatile distributions and lower for stable plumes).The exponential energy function E associated with the Gaussian factor f o for a single z k,i is therefore given as , where ∆t z k,i is the time difference between the current time step and the time the measurement z k,i was recorded.The factor for each node x i with a nonempty set of measurements Z i can be expressed as where |Z i | refers to the cardinality of the measurement set associated with the i-th node and ).The regularisation factor f r defines the potential between the connected nodes in the graph and relates to the influence that the gas concentration at voxel i has on each its neighbours j ∈ N i , given as l i,j = x i − x j .In the 3D case, the regularisation factor is also defined between voxels that are above and below x i (as shown in Fig. 1) and thus |N i | ≤ 6.The associated Gaussian for regularisation is modelled by l i,j ∼ N (0, σ 2 r ).By tuning σ 2 r relative the σ 2 o , we scale how much a measurement will affect its surrounding nodes.Regularisation factors also account for the effect that an obstacle between nodes has on the distribution by conditioning the factor graph based on an occupancy map.Let o ij ∈ {0, 1} refer to the binary occupancy between node i and its neighbour j.When o ij = 1, there is an obstacle blocking the correlation of the two cells' concentrations.If o ij = 0, free space exists between the nodes, and therefore the regularisation factor acts uninterrupted.This gives the governing energy equation for regularisation of a pair of nodes (l i,j ) that are conditioned on the occupancy map as E r (l i,j ) = 0.5( ).Thus, the regularisation factor associated with node i can be written as where ).Note that the regularisation factor contribute to both the edge potentials and the selfpotential of the related nodes.
Finally, the default factor f d is defined which lightly anchors the concentrations at cells to a nominal level (z 0,i ), given by the Gaussian . This is required to enforce the cell's concentration close to zero (or any other arbitrary background concentration value) in the absence of any information in its surrounding area.If this is not set, then sensor measurements can propagate infinitely far into the factor graph causing high predictions in unsampled regions.This factor will also improve the numerical stability for isolated local cells.By tuning this factor relative to f o , we can scale how much a sensor measurement overrides the zero prior value.This factor is similar to f o in that it contributes to a cells self potential.The energy related to the default factor for node x i is given as E d (z 0,i ) = 0.5(x i − z 0 ) 2 /σ 2 d , and thus the default factor for all nodes can be expressed as Based on the established definitions of the individual factors, the information matrix Λ and the observation vector g of the GDM problem can be constructed.Collecting the first term in ( 4)-( 6), it can be derived that the diagonal entries of the information matrix Λ follows whereas the second term in ( 4) and ( 6) contributes to the i-th entry of g such that The off-diagonal entries of the information matrix Λ ij depend on the connectivity of the graph.The non-zero entries are associated to the edges of the graph and can be formed as As ( 7) -( 8) fully define the 3D-GDM GMRF model, an appropriate solver is then needed to provide the solution to the inference problem.The marginal inference problem for each cell focuses on p(x i ) = N (µ i , P −1 i ), where µ i = {Λ −1 g} i is the mean and P i is the inverse of the variance (i.e., precision) of node i, given by P It should be noted that with the increased size of N i to account for the 3D problem, inference will be performed on an information matrix significantly less sparse than the 2D case.Furthermore, unless a larger discretisation size is chosen, |X | will also significantly increase.Therefore, the combination of these two characteristics means that solutions with direct solvers very quickly become infeasible.For large scale environmental monitoring, more efficient solvers must be implemented in order to retain online inference despite the scale of the problem.
In this implementation we model gas distribution characteristics using only concentration measurements, despite both [5] and [8] showing that including anemometer measurements can provide more accurate maps.Our main focus experimentally is on indoor locations where the wind field is not significant enough to dominate gas transport, so we instead optimise the solver for GDM, leaving the extension of modelling the wind field as future work.

III. GAUSSIAN BELIEF PROPAGATION FOR GDM
In [7], after establishing the GMRF model, the solution of the inference is then calculated as the convenient least squares problem which in the original work is solved using the Cholesky decomposition to solve the Gauss-Newton method equations.There are many different solutions to inference on a GMRF, as discussed in [15], and by utilising a solver that suits the application of a mobile robotic sampling agent, real time estimates of complex distributions can be achieved.
Belief propagation is a general tool for inference on a factor graph and functions by passing iteratively updating belief messages between neighbourhoods in the factor graph.This message passing process is equivalent to Pearl's local messsage-passing algorithm [21].Due to the Gaussian assumption placed upon the factor graph, a message passed between nodes is a Gaussian belief.The following description is based the procedure by Shental et.al [22] and an excellent visual overview of belief propagation by Ortiz et.al can be found in [23].A message consists of real values that are passed between nodes and are governed by two rules: the 'sumproduct rule' and the 'product rule'.Since we are dealing with a Markov random field of continuous variables, the sumproduct rule becomes the integral-product rule [24].Therefore, for a message sent from node i to node j the following applies: (10) where the notation N i \ j is the set of neighbours of i not including node j.The marginals are computed as per the product rule: where γ is a normalisation constant.Fig. 2 shows an example message passing between node i and node j in a small neighbourhood, with symbols corresponding to equations (10) and (11).The self potentials established in section II-A are of Gaussian form so that φ i ∝ N (µ ii , , P −1 ii ) and m ki ∝ N (µ ki , P −1 ki ).Then based on the rules of the product of Gaussians shown in [22], the Gaussian associated with right hand side of equation ( 10) can be calculated as: where P ii Λ ii (a-priori precision) and µ ii g i /Λ ii (mean of the self potential).After resolving the remaining Gaussian integral of equation (10), the message m ij (x i ) that is sent along an edge for any given pair of nodes in the network is given as Eq. ( 14) shows how updated beliefs are sent from one node in the network to a neighbour, and then the next step is to repeatedly send messages between nodes under a specific scheduling until the messages reach a convergence threshold.The final step is to calculate the marginals of all x i ∝ N (µ i , P −1 i ) ∈ X as per equation (11) to give: This completes the GaBP process and shows that since messages are passed locally with no specific ordering and that marginal mean and uncertainty values are calculated iteratively, the solver can flexibly resolve in any local region of the graph and recover the current solution vector at any time.This is a key motivation for using the GaBP solver in the mobile robotic domain wherein information is usually acquired within a local region and can be without a fixed frequency.
Comparing the estimation accuracy of such iterative solvers against their direct counterparts, the findings of Weiss et al. [24] show that the GaBP solver will converge to the exact marginal means estimated by direct methods, even when the graph exhibits many loops (which is strongly the case in GDM given the interconnected nature of the spatial distribution).However, in these very loopy graphs the variances will often be incorrect (always overconfident in GMRFs), but as stated in [24], for the marginals to converge to the exact values then the variances must be heuristically correct, indicating that despite their real values being incorrect, the relative values of the variances between nodes in the network must be correct.
To further increase computational efficiency, the number of messages sent by GaBP can be reduced from O(n 2 ) to O(n) per iteration round when broadcasting aggregated sum of potentials for each node, as per Algorithm 1 of [22].Since GaBP is a general tool for solving factor graphs, the fact we are applying the solver to a GDM problem or that it is a 3D-GDM problem has no effect on its fundamental application.Therefore, equations ( 12) -( 16) can be directly applied to the 3D graph {Λ, g} derived in section II-B.

IV. 3D-GDM SOLVER
Based on the formulation of the 3D factor graph and the principles of general Gaussian belief propagation, the proposed algorithm to achieve fast 3D inference will be outlined in this section.Key improvements over the current state-of-the-art include a hybrid message scheduler as well as a dynamically built factor graph based on an information theoretic criterion.With these changes in place, efficient local 3D-GDM can be achieved.

A. Hybrid schedule with sporadic measurements
In the introductory work on G-GaBP [16], several message scheduling algorithms are tested in a GDM scenario.These are random, round-robin, residual [25] and wildfire [18] schedulers.It is found, for the specificity of environmental monitoring with a point sensor, that the wildfire algorithm has much better mean and variance convergence around the measurement region than the other schedules.It has also been found however that in unsampled regions away from the local vicinity, that the speed of uncertainty convergence is much slower when using the wildfire algorithm.In contrast, residual belief propagation shows good uncertainty convergence globally but with a poor rate of marginal mean convergence.Furthermore, when performing wildfire message propagation, messages are sent until a threshold information theoretic criteria is met, which leads to variable iteration times (depending on the amount of new information a particular new measurement contributes to the factor graph).Thus, there may be unused computation time between sensor readings where the inference has adequately converged in a local vicinity.This is especially relevant in mobile GDM where concentration readings are often subject to a transient response [26] and agents may need to stop and sample a location for a period of time [27], [28].It would be useful in this downtime to continue resolving the factor graph in order to maximise the online performance of the mapping agent.Based on these insights, we propose a new hybrid schedule that leverages the rapid local convergence of the wildfire algorithm, combined with the global convergence of the residual belief propagation schedule.
To achieve this, we enforce a mode switching that changes the message schedule to residual propagation when a wildfire iteration finishes.Since both wildfire and residual propagation track the same message residuals, extra computation is not required and the solver can simply switch between modes given its current convergence.The residual of a message is defined as the Bhattacharyya distance between the current message being sent and the previous message sent along the same edge: Other information theoretic metrics can also be used and as discussed in [16], but performance of the solver is more affected by the schedule rather than the specific metric (for the case of GDM).Formally, the hybrid message scheduler is outlined in Alg. 1.
Algorithm 1: Hybrid message schedule Update and send m tj according to Eq. ( 12)-( 14) ∀j ∈ N t Update and send m tj according to Eq. ( 12)-( 14) Due to the nature of belief propagation wherein the marginal mean and precision values are calculated as part of the message passing process, the current global marginal means and variances can be extracted at any time during both the wildfire and residual loops.The efficiency gain of using a hybrid message schedule is tested in simulation in section V.As stated in Introduction, this addition is a key contribution to efficient GDM and is combined with our third contribution to the algorithm, to be described in the following subsection.

B. Local graph building for local inference
Thus far, our work on G-GaBP has focused on how to utilise an efficient local solver in order to overcome the computational inefficiencies of the direct methods.Despite the shift to local inference, until now the factor graph upon which we perform inference is still in a global frame and considers the entire area to be searched.This is naturally opposed to the local nature of GaBP and is something we wish to address.
In the works of [7], [8], [16], [29], the area is first discretised into cells and then the factor graph is built in its entirety upon this a-priori discretisation.This means that even with GaBP, inference is being performed in areas far away from the sensor location and efforts are wasted in calculating these distant marginals which approach {µ i = 0, P i = P ii } i.e. the a-prior state.Furthermore, significant computation is involved with keeping the connections of such a large factor graph up to date i.e. graph management (based on updating occupancy values).Therefore, it is pertinent to consider a local graph construction scheme that iteratively redefines the inference problem dependent on the expected amount of information gained by performing said inference.This gives rise to a chicken and egg problem, where we do not know the problem formulation before performing inference, something which is impossible with conventional direct solvers.However, the message passing of GaBP does allow us to achieve this.The rest of this section will outline our third contribution, which is a graph management methodology for redefining a local factor graph whilst inference is being performed upon it.
To achieve iterative graph construction, we can take inspiration from iterative factor graph problems such as those generated in pose optimisation for SLAM.In these works at each iterative step, new sensor measurements and landmark observations are added sequentially to the factor graph and connected to a new corresponding state estimate and/or a previous state estimate in the case of loop closure.Whilst this process can be mirrored in GDM by adding a new graph node for each new sensor reading at a previously unvisited location, it is not immediately obvious how far this sensor measurement will effect is local surroundings.It would be feasible to heuristically add a graph support to the initial voxel location, but this would not take into account the current belief and therefore the information actually added by the said measurement.We can however take advantage of the properties of wildfire propagation to expand our factor graph in an information theoretic manner, as messages propagate out from the measurement.The following outlines the process of this graph expansion.
When a new measurement z k,i is recorded from a sensor, whose corresponding node x i / ∈ X , a new variable node is added to the factor graph i.e.X + ← − x i , registered at the given sensor location.The six possible neighbours of x i are then added to the neighbourhood of x i given there is no obstacle blocking the connection between any pair of voxels.The information matrix Λ and observation vector g are updated using Eq. ( 7) -( 8) to account for the new nodes in the factor graph.The wildfire iteration can then occur with some modification to allow for dynamically expanding the factor graph around the new information.As per the standard wildfire process, when a message is sent to neighbour j from node i, if the residual is greater than the wildfire threshold r m (m ij ) > , the neighbour j is added to the queue Q + ← − j for subsequent transmission (Alg. 1, Lines 8-10).However, at this point the neighbour j may not have been expanded i.e., N j = i, and therefore cannot propagate any new information further.In this case, the node j should then be checked for potential expandable neighbour nodes, e, to add to the factor graph X + ← − x e , ∀ e ∈ N j ∧ x e / ∈ X .This process continues as the wildfire propagates and the factor graph grows iteratively until the outer nodes of the factor graph are receiving less information than the wildfire threshold .This then couples the size the graph to the information theoretic value of a given measurement.Due to the fact that both information propagation and graph construction depend on the value of , selection of this value is important for efficient convergence and adequate graph coverage.
Fig. 3 shows an example 2D iterative growth process pictorially.In this scenario the measurement z k is added to a discrete point in 2D space and an associated node x i is generated along with 3 neighbours {x j1 , x j2 , x j3 }.Note that a possible neighbour to the left has not been added due to the presence of an obstacle between x i and the potential x j4 location.The information from x i is then sent to each of the neighbours, wherein both x j1 and x j2 receive information above threshold .These nodes are then added to the queue Q + ← − {x j1 , x j2 }, and then their respective neighbours x e ∈ N j are added to the graph (green nodes) if x e / ∈ X .Conversely, x j3 does not receive above the threshold information gain, and therefore its potential neighbours are not added to the graph.Also note, that a independent part of the factor graph (purple nodes) from a previous measurement (z k−1 ) may also be connected to the new expansion and will iterate messages as per Alg. 1.If the message information does not propagate enough to connect to the main graph, then the independent sub-graph may also remain disconnected forming two independent graphs.Due to message passing only needing neighbour information, having disconnected graphs makes no difference to the global solution and can be reconnected at any time.
In the circumstance where a new measurement does have an associated node z k,i , then belief propagation occurs as normal until an unexpanded node is met, after which its neighbours are attempted to be added as per above.

C. Proposed algorithm
As the proposed methodology shifting away from an apriori graph definition towards a dynamic inference problem, we now introduce the proposed G-GaBP algorithm, whose pseudo-code is provided in Alg. 2.
Firstly, hyper-parameters of the problem are set a-priori and the empty factor graph is initialised.When a new measurement is received, if z k,i has no corresponding x i ∈ X , then a new node is added along with its immediate neighbours (lines 6-8).When adding these connections to the factor graph, new edge potentials are initialised which define the incoming information between the new nodes (line 9).During this step, the prior message information, m k−1 ij , must be set to some minimum amount of information that can be received along an edge.This is defined as N (0, σ 2 p ) by studying Eq. ( 14) and the case of a factor graph consisting of a chain of nodes at an infinite distance from an observation.This is required since if this value is initially set to {0,0}, then calculating Eq. ( 17) results in infinite information and therefore the propagation of growing the factor graph will occur indefinitely until the entire area is represented in the factor graph (nullifying the local graph growing procedure).
The information matrix Λ and the observation vector g are then updated to reflect any new nodes added to the factor graph and also to recalculate the time dependent precisions of Eq. 7 and time dependent observations of (8) (line 11).The self potentials required for belief propagation are then also updated to reflect this change (line 12).
Line 13 then starts the wildfire propagation process.The node associated with the latest measurement i.e., z k,i , is first added to the queue for propagation.Propagation then occurs according to the wildfire schedule, with the factor graph iteratively expanding with the propagating information (lines [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32].If the time dependent decay factor has been employed, after the first measurement has propagated its information, then all historic measurements must also trigger their own wildfire iterations to transmit the now reduced information which they carry.Since residuals have already converged around previous measurements, these subsequent wildfire iterations usually require very little message passing.
Once all measurements have converged, residual belief propagation is performed until a new measurement is received Add new edge potentials: Update Λ and g using ( 7)-( 8) Update self potentials φ i : Update and send m tj according to 18 Eq. ( 12)-( 14) Update and send m tj based on Eq. ( 12)-( 14) 37 end (lines 35-36).If a new measurement is received before this stage, the process can be stopped at any point and the new measurement inserted (the anytime property of GaBP).This concludes the algorithmic portion of the paper.The proposed solver of Alg. 2 has also been developed in python as an open source ROS package 1 for implementation on real and simulated mobile robots (for both 2D and 3D inference).

V. SIMULATION STUDY
To showcase the real-time application of the 3D-GDM solver for use with a mobile robot, we deploy the algorithm in an example online mapping scenario.The setting for the simulation is a time-varying source release in a 200 × 100 × 20m industrial arena consisting of buildings and shipping containers.The ground truth concentration data is modelled using GADEN [30], a 3D filament dispersion simulator capable of modelling gas dispersions in feature rich environments.The specific release for the scenario is an acetone source situated in CFD derived wind flow field with release rate of 1000ppm/s.Whilst the source release is modelled as an unsteady process, due to the constant wind field, the main structures of the plume remain consistent and therefore the mapping process is approximately time invariant.The sensor response behaviour is modelled inside GADEN and is inherently noisy (see [30] for further details on the simulator).The simulated mobile sensing agent is a quadrotor with a photoionisation detector (PID) affixed to the base.Realistic physics and flight control are handled using the popular simulation environment Gazebo and the PX4 autopilot plugin respectively.The operating scenario of the mapping task is given in Fig. 4.

A. 2D simulation study
To benchmark the algorithm against the direct methods, we first deploy the algorithm in a 2D inference task over the entire area.This is necessary due to the scale of the area making direct solvers incapable of performing 3D inference.For example, at 1m resolution the size of the corresponding GMRF for the 3D simulation scenario is a 4 × 10 5 -by-4 × 10 5 information matrix, which takes approximately 6 minutes to resolve using a normal PC.This is clearly unsuitable for online mapping tasks.We also make comparison to our original G-GaBP algorithm [16], without the optimisations proposed for the 3D-GDM solver.To differentiate between those two 1 GaBP open source ROS package available from https://github.com/callum-rhodes/gabp mapping algorithms, we dub the proposed solver as GaBP+ and the original as GaBP.
For data collection, the unmanned aerial vehicle (UAV) flies at 1m/s along a pre-determined sweeping trajectory through the plume area with fixed sampling waypoints that must be sampled (at 3m separation distance).If the UAV encounters an obstacle, it performs obstacle avoidance manoeuvres and any measurements taken above the estimation plane whilst avoiding obstacles are not included.The concentration sensor publishes at a constant rate of 2Hz and the mapping algorithm can include measurements that are captured when travelling between the waypoints.It is also enforced that at a minimum, the UAV must resolve the measurement at each fixed sampling point into the distribution map, before sampling the next location.This constraint is chosen to represent the scenario where a decision about the next course of action may be taken based on the updated belief of the dispersion (as would be required by informative path planning or dictated by a human operator).It also allows the dataset Z to cover the search area evenly.This arrangement means that once a sample is taken at a waypoint, the measurement is inserted into the solver and whilst resolving said measurement, the UAV will start to move to the next waypoint.If by the time the UAV reaches the next waypoint, that the previous sample has not been resolved into the map, the UAV will wait before sampling the current location.Conversely, if the solver is capable of resolving the first measurement before reaching the next waypoint, then any intermediary measurements taken en-route may be included in the inference as soon as the solver becomes available to accept a new measurement, i.e., sensor measurements must be sequential.For GMRF this is defined when the MAP solution of p(x|Z) is retrieved, and for GaBP this is defined as when the wildfire iteration for all current z k,i ∈ Z completes.The ability to include intermediary measurements is chosen to demonstrate the effective online performance of the GaBP based solvers.The ground truth data is also gathered for error comparison, which records the average concentration over the simulation at each 1m grid cell over time.For context, if a single mobile sensor was to collect the ground truth data with a 10s sampling period, the data collection would take over 12 hours for the 2D ground truth (and over 60 hours for 3D), whereas the expected search time for the 2D sweep is approximately 1700s (within typical flight budgets for quadrotor UAVs).From an application point of view, this shows the necessity for such dense mapping algorithms to exist for hazardous material first response.
Hyper-parameters for the GMRF model are kept consistent between the 3 algorithms with σ 2 r = 2, σ 2 s = 0.1, σ 2 d = 1e 4 and for the GaBP solvers, = 0.01.The time varying factor is set as σ 2 ζ = ∞ since we are modelling a quasi-steady state plume structure and therefore measurements are time independent (but still subject to noise).Fig. 5 shows the resultant marginal mean concentration maps after completing the search, and Fig. 6 the overall reduction in root mean square error (RMSE) over time, for each of the 3 algorithms.Numerical statistics of the mapping task for each algorithm are shown in Tab.I. Metrics for comparing uncertainty outputs between GMRF and GaBP are not included since their output values are not comparable, as shown in [16].
Only the region of the interest (the plume) is used for the RMSE calculation, defined as the set of cells, X plume , exhibiting higher than a background level concentration, z thresh = 100ppb, i.e., X plume ← x i ∈ X | µ gt (x i ) > z thresh where µ gt (x i ) is the time averaged ground truth concentration at the location associated with variable x i .This is so that the noise associated with modelling the zero concentration background regions does not dominate the output RMSE of each algorithm.The calculation for RMSE is therefore: where n plume refers to the size of the set X plume .The resultant marginal mean concentration maps show very similar structures of the plume to the ground truth, with variations in peak values (particularly near the source) due to the chaotic nature of gas dispersion and a short sampling time.Whilst the mean maps show similar features, when appreciating the RMSE relative to the ground truth it can Fig. 6. 2D simulation RMSE reduction over time for each of the 3 tested algorithms be seen that G-GaBP+ has the best final state estimate and error reduction over time due to the hybrid solver taking advantage of downtime between wildfire iterations.Because of this, we also see that measurements are also resolved slightly quicker (14ms vs 18ms), since they are being inserted into a more converged belief.This improved performance is managed whilst only estimating one quarter the number of states (6536 vs 25228), significantly reducing the computational overhead in managing the factor graph structure, which is linearly related to the number of estimated states.The direct GMRF solver in contrast to the GaBP solvers has much worse online performance.In our simulation, this is directly caused by the slow update rate of the GMRF solver and therefore significantly fewer measurements are recorded in the final map (460 vs 3253).Whilst GMRF takes an average of 4.8s to resolve the map per measurement, GaBP+ is capable of solving a single measurement on average in 14ms, significantly quicker than the update rate of the sensor and hence why GaBP and GaBP+ process the same number of measurements (i.e.all of the measurements output by the sensor).Reaffirming our conclusions in [16], GaBP propagation has shown its superior performance for online inference against direct solvers.The following simulations explore the online obstacle aware 3D mapping abilities of the proposed algorithm, something not yet achieved in the literature.

B. 3D simulation study
For the 3D study, we explore the same scenario however expand the state estimation to include the 3D advancements of section II-B.This expands the possible estimation set X , to be up to 4×10 5 nodes which are also more densely connected than the 2D case.As explained previously, GMRF is unsuitable for online state estimation in this scenario and therefore GaBP and GaBP+ are compared together.The graph management overhead for the original GaBP algorithm is also computationally limited when considering the full 200×100×20m scenario area and therefore only the area covered by the sweep is considered in the 3D task (furthermore showing the necessity of the contribution in section IV-B).Parameters of the search from the 2D study are carried forward however a saw-tooth sweep in the vertical axis z= {1, 3, 5}m is performed to collect measurements required for 3D mapping (see Fig. 4).This extends the expected search time to 70 minutes.The mapping results are given in Fig. 7, which shows 1m slices of the estimated 3D map compared against the ground truth data.Fig. 8 shows the RMSE reduction over time.Other Statistical data are provided in Tab.I.
As with the 2D task, concentration estimates across the map are visually comparable with the ground truth and the difference between GaBP and GaBP+ is minimal.However, when viewing the RMSE reduction rate it is clearly shown that the proposed 3D-GDM algorithm shows a numerically greater estimation accuracy over time, and does so with far fewer number of estimated states.As with the 2D case, the average resolve time per measurement is very low for both algorithms (less than the 2D case) with GaBP+ being marginally faster than the original as previously discussed.This is because the plume above 3m is not very different from the a-prior state (µ = 0) and therefore measurements taken above 3m are close to 0 and do not take long to converge (reducing the average measurement resolve time).This is somewhat counter-intuitive as it would be expected that a 3D mapping task would be more complicated and take longer to converge.However, wildfire propagation is scale agnostic and its information theoretic nature allows areas of low information to not negatively affect the solve time.Furthermore, the impact of obstacles on the 3D gas distribution is fully modelled as zero concentration readings are predicted in the obstacle locations at heights of z= {1, 2, 3} but then smooth predictions are shown in these same areas for variables at height z= {4, 5}.
This concludes the simulation study and it has been clearly demonstrated the advantage of the GaBP over conventional direct solvers.In addition, the locally expanding inference and hybrid schedule proposed in this paper have demonstrated a significant improvement on estimate convergence as well as on the computational efficiency of graph management.These additions have decoupled the problem definition from a predefined search area allowing for iterative, fast and local 3D distribution mapping, a major step forward for complex online GDM for mobile sensing.

VI. EXPERIMENTAL TRIAL
Having shown the benefit of the proposed 3D-GDM solver against the current state of the art in simulation, we now apply the algorithm in a real-world 3D mapping task to show its capability.One of the key benefits of using the GaBP solver is the ability to generate real time mapping solutions onboard the local hardware of the sensing agent, whereas all work previously has been performed offboard.This section will demonstrate online 3D mapping entirely onboard a ground vehicle (with networking only used for monitoring), and is the first practical application of the GaBP solver for GDM.The scenario for the task is a cluttered industrial workshop with several isolated work units within the area.Three source releases are placed within the environment at different heights (this is to discern the concentration levels along the z-axis).The source releases are beakers of liquid acetone, agitated by a ultrasonic bubbler and dispersed using a desk fan.
To perform data collection, a Clearpath Husky unmanned ground vehicle (UGV) is fitted with a sensing pole that has PID sensors affixed at 0.3, 1 and 1.8m.Sensor data is recorded at 1Hz for each sensor and fed to the GaBP solver in an adhoc manner.The UGV is manually teleoperated around the environment in an exploratory search with live feedback of the marginal mean and uncertainty maps.Online 3D obstacle maps are built using the Lidar outputs of a Velodyne VLP-32C and a Robosense-Bpearl to achieve nearly full 360 o coverage of the surroundings.The 3D map is generated using the popular 3D map representation and ROS package Octomap [31], and it is occupancy values from this software that directly condition the factor graph structure defined in Eq. ( 7) - (8).Localisation of the robot is achieved using the Velodyne VLP-32C in conjunction with the open source graph SLAM package hdlgraph-slam [32].Good localisation is fundamental for accurate GDM mapping, not only because pointclouds are inserted into the Octomap relative to the local pose (i.e.mapping with known poses), but also the pose estimate informs the measurement factor association with a specific state variable (i.e.defines z k,i ).The 3D G-GaBP algorithm is run with the same hyper-parameters as in the simulations however with a higher resolution of 0.5m and a higher = 1.The approximate size of the area is 40 × 25 × 10m.
Contrary to many other GDM works, the proposed experimental platform does not rely on a external positioning system or an offboard processor, and therefore the platform can be deployed in any scenario in complete isolation.To the best of the authors' knowledge, this is the first time such a versatile platform has been deployed for the task of GDM, aside from the academic contributions of the proposed algorithm.For a video of the 3D-GDM experimental trial, please visit https://youtu.be/crAJd4afW8c.Fig. 9, shows the 3D plume predicted by the proposed algorithm for each one of the 3 sources as the robot is driven past each source.The 3D aspect of each sources plume can be easily discerned when studying the coloured voxel concentration map.Source 1 is the first source encountered and is situated above all three sensors at 2.3m.A high concentration is measured next to the source location (bright yellow voxel) but nearly zero chemical is detected on the bottom sensor at the same location, this results in a plume that is dispersed above the floor settling in the opposite room at floor level.We can also see significant mixing behind the source at the back of the room.Source 2 is encountered next, and this source is dispersed along floor level at 0.2m.Only the low sensor reads a high concentration value and therefore the distribution shows the plume spreading along the floor, with some mixing at 1m.To reinforce the importance of the 3D modelling aspect, a 2.5D map is also shown in the top left corner at the point in time the source is sampled.To generate the 2.5D map, concentration values from the 3D inference are averaged along the z-axis and projected down to their corresponding x-y coordinate.If using this map alone to analyse the gas distribution, then the high altitude source dispersion is clear, however the low height source is masked by the low concentration readings of the mid and high sensors and therefore cannot be discerned as a possible high concentration area (whereas this source is clear using the 3D plume map).This is a significant advantage of using a 3D inference model over a 2D model since source 2 may have been overlooked if using the 2D map alone.Finally, source 3 is sampled in a larger room and the plume is shown to disperse up to about 1.5m high.The low concentration readings of the high sensor stop the distribution propagating upwards where the acetone does not mix (due to a higher density than air).The complete distribution map is shown at the bottom of Fig. 9 after the robot has searched the area.Due to the density of vapourised acetone and its propensity to settle on surfaces, we can see that the low and mid sources do not show as much dispersion as the high source, which is able to be transported much further through the environment.From this image we can also see the requirement for obstacle aware GDM in such a cluttered environment, since for both source 1 and source 3 there is no mixing of the high concentration to the other side of their surrounding walls (and equally the low concentration readings outside the rooms do not directly interact with the high concentration plume inside the rooms).
Mapping statistics for the entire experimental mapping task are shown in Tab.II.Due to the fact that several returning passes of the rooms are made during the 15 minute search (and a higher value), the average resolve time per measurement is significantly lower than the simulation tests since the main plume structures are established within the first pass, and then only slightly updated and refined in subsequent traversals.

VII. FUTURE WORKS
The work encompassed in this paper has demonstrated the real-world applicability of the belief propagation solver to an environmental monitoring task with a mobile sensor.Having validated that the solver can perform large scale, real time and obstacle aware modelling onboard a standalone platform, efforts can now be made to further develop the system and make use of the computational advantages brought forward by GaBP.One of the key areas this work can be extended is in more complex modelling tasks.A natural extension is to included wind measurements alongside concentration measurements, and to model both phenomena and their mutual influence concurrently, such as that proposed by Gongora et.al [8].It is foreseen that belief propagation will excel in this case, as GaBP has already shown to provide excellent solve times for other non-linear optimisation tasks, as shown by Ortiz et.al [23].Furthermore, predictive inference for informative path planning (IPP) using the GaBP framework now becomes computationally achievable, and with the ability to have a fixed time cost.However, due to the the convergence of estimates also factoring into the quantified uncertainty output, how to apply information theoretic criteria remains an open question.Finally, distributed systems may be explored wherein a multirobot system uses the GaBP framework to update neighbouring robots of shared beliefs in an efficient manner.

VIII. CONCLUSIONS
This paper first outlines the problem of performing inference for 3D gas distribution maps upon a connected factor graph.Due to the computational complexity of such a large problem, the GaBP algorithm is applied as the iterative solver to achieve real time distribution maps as concentration readings are collected.The G-GaBP algorithm is then extended over the original of [16] in two key areas.First, a hybrid message schedule is proposed that leverages the fast local convergence of the wildfire algorithm with the globally converging properties of residual belief propagation.Second, the GDM problem is formulated incrementally as inference is performed, tightly coupling information propagation and graph construction in order to achieve entirely local inference.This feature also allows the proposed algorithm to be more easily integrated with popular SLAM solutions.
These new contributions are then compared in simulation against the standard direct solver of G-GMRF and the original G-GaBP in both 2D and 3D mapping tasks.The proposed solver is found to be able to perform more accurate inference (due to the hybrid scheduler) whilst simultaneously estimating significantly less states (reducing computational burden on factor graph management).These enhancements allow the 3D G-GaBP algorithm to be deployed onboard cheaper hardware units without the need for GPU acceleration.
To prove this and show its real-worlds applicability, the algorithm is deployed onboard a UGV with 3 gas sensors at varying heights.The vehicle is driven around a cluttered GPS-denied environment and in real time resolves complex structures of the gas distribution for 3 different sources placed at different heights in the environment.It is also shown how using a 2D distribution map can obfuscate features of the plume that are only identifiable using 3D inference.The development of an online 3D-GDM solution opens up further

Fig. 1 .
Fig. 1.Factor graph cliques for each type of factor in the 3D-GDM problem.

Fig. 2 .
Fig. 2. Local message passing from node i to its neighbour node j.

Fig. 3 .
Fig. 3. Information theoretic local factor graph growth around sequential measurements

Algorithm 2 : 4 if 5 if x i / ∈ X then 6 X+← − x i 7 Define N i based on o ij 8 X+←
Proposed G-GaBP solver 1 Define factors f p , f o , f d and z 0 = 0 2 Initialise empty factor graph, X , Λ, g ← ∅ 3 while Active do New z k,i then // Update factor graph − x j , ∀j ∈ N i 9

Fig. 4 .
Fig. 4. Simulation scenario.2D sweep is shown left with black trajectory and fixed sampling waypoints (circles).The green circle and arrow shows the source location and wind direction respectively.Pictured right is the simulated UAV operating with a representative 3D sawtooth sweep (red circles and trajectory) for the 3D mapping task.Also shown is the online Octomap being built as well as instantaneous gas particles from GADEN (green dots).

Fig. 5 .
Fig. 5. Final marginal mean estimates of the plume region for each algorithm compared against ground truth values.Obstacle locations shown as grey rectangles.

Fig. 8 .
Fig. 8. 3D simulation RMSE reduction over time for GaBP+ and GaBP.Also shown is the number of estimated variables, |X |, over time.

Fig. 9 .
Fig. 9. Experimental application of the 3D G-GaBP solver in a cluttered GPS-denied scenario.Each estimated 3D plume structure for the 3 sources are shown (coloured voxels for any cell with a concentration above 20% of the maximum concentration value, yellow indicates high concentration).The corresponding 2.5D map at the point where the source is driven past is shown in the top left corner.The base of each red arrow indicates the source location and the direction represents the fan induced wind direction.The real-world location of each source plume is shown on the right for context.Shown bottom is the full 3D & 2.5D distribution maps superimposed on the Octomap.

TABLE I MAPPING
STATISTICS FOR EACH ALGORITHM FOR BOTH THE 2D AND 3D SIMULATIONS