Determining the three-dimensional structure of a volcanic plume using Unoccupied Aerial System (UAS) imagery

a Department of Aerospace Engineering, University of Bristol, Bristol BS8 1TR, UK b Department of Computer Science, University of Bristol, Bristol BS8 1UB, UK c School of Earth Sciences, University of Bristol, Bristol BS8 1RJ, UK d Department of Earth Sciences, University of Cambridge, CB2 3EQ, UK e Instituto Nacional de Sismología, Vulcanología, Meteorología e Hidrología (INSIVUMEH), 7a Av. 14-57, Zona 13, Guatemala


Motivation
Volcanoes release gas, ash and aerosols into the atmosphere, in the form of either passive or explosive emissions. The ash-rich plumes from explosive eruptions can reach the altitudes of commercial aircraft flight paths. Here, ash particles pose a hazard to both the airframe and engine, and thus require avoidance strategies that can generate regional and global socio-economic impacts (Chen and Zhao, 2015) (Guffanti and Tupper, 2015). Airline authorities currently have insufficient information to confidently ensure the safety of jet aircraft (Prata, 2009) and require more detailed knowledge of the 3D distribution of ash in near real-time to re-route aircraft efficiently.
Currently, volcanic emissions are mainly monitored using a combination of satellite-based observations (Thomas and Watson, 2010) and dispersion modelling (Peterson et al., 2015). Dispersion models can provide 3D forecasts of volcanic cloud trajectories (Stohl et al., 2011) that are 'fully volumetric' (i.e. the inner regions of plumes are modelled) and up to hemispheric in scale. However, these model algorithms are initialized based on plume source parameters (Mastin et al., 2009). Poor constraints on these initial parameters can therefore result in large uncertainties in the spatial distribution of volcanic plumes and their physical properties (Devenish et al., 2012) (Mastin et al., 2009) (Beckett et al., 2014). There is a need, therefore, to collect plume source parameters (which may be at the meter scale) in order to initialize the plume dispersion models (which may be km in scale).
Advances in Unoccupied Aerial Systems (UAS), combined with three-dimensional computer vision methods, provide interesting measurement opportunities for monitoring initial plume properties. It should be noted, that the magnitude of the volcanic event to be observed will be somewhat proportional to the size of UAS required to observe it. Commonly deployed UAS, similar to that used in this study, will only be suited to mapping of decimetre-scale plumes; although the method could be applied to imagery from larger UAS, or even manned aircraft, as required. Computer vision methods have already been applied to 3D mapping of volcanic plumes from the ground using space carving methods , but this previous work used ground-based cameras and a non-probabilistic approach to the space carving. Space carving has also been applied to satellite images using Structure from Motion (SfM) methods (Zakšek et al., 2018). The space carving technique used in these articles does not produce fully 3D volumetric models. Instead, it generates a 3D external surface or 'visual hull' using an averaged position for the outer optically thick depth. This is what is meant by the term '3D reconstruction' used throughout this paper. Cameras mounted on UAS offer an alternative viewing geometry to ground observations and at a higher spatial resolution than satellite imagery (Remondino et al., 2011) (Colomina and Molina, 2014) (Harvey et al., 2016) (Beni et al., 2019) (Gomez and Kennedy, 2018) Some initial UAS-based 3D measurements of volcanic plumes have already been achieved (Gomez and Kennedy, 2018).
The method presented in this study has been developed using a case-study dataset acquired at Volcán Pacaya, Guatemala. However, in principle, the process can be applied to other plumes and UAS camera systems. The additional insight into plume topology provided by this method bridges the gap between satellite and ground-based monitoring techniques. The images capture a sub-section of the time evolving condensed steam plume emitted from the main vent of Pacaya. Since the entire plume is not captured in the images, the purpose of this study is not to provide quantitative measures of volcanic behaviour, but rather to demonstrate proof of concept. The ultimate motivation of the current work is to provide a method for the measurement of volcanic gas/ash plume 3D spatial properties using UAS-mounted multispectral thermal imaging.

Volcanic plume measurements
Ground-based measurements of volcanic cloud properties can be made by weather radar (Lacasse et al., 2004), Doppler radar (Hort and Scharff, 2016), LiDAR (Hervo et al., 2012), or visual/multispectral cameras and spectrometers, e.g. (Lopez et al., 2015). Due to the oftenchallenging conditions (danger and inaccessibility) around active volcanic peaks, such sensors are normally stationary and positioned far away, often only providing views from below or level with the plume. Satellite imagery is used extensively to monitor volcanoes and mitigate risks of their hazards as it provides wider coverage of the Earth (Zehner, 2012). Several satellite instruments have proven their capabilities in detecting most hazardous volcanic clouds (Thomas and Watson, 2010). These instruments can provide near-global coverage at all times, but their main limitations are temporal and spatial resolutiontypically providing a single synoptic snapshot. UAS are now starting to find utility in volcanology where they are being used for the in-situ measurement of volcanic emissions (McGonigle et al., 2008) (Pieri et al., 2013) (Stix et al., 2018) (Shinohara, 2013) (Mori et al., 2016) (Liu et al., 2019) topographic mapping (Rokhmana and Andaru, 2016) (Darmawan et al., 2017) (Carr et al., 2018), and identification of volcanic instabilities (Darmawan et al., 2018). UAS observations are also not without their challenges, however. Even with high acquisition frequency, it takes time to make enough observations from various positions to undertake a reconstruction. With a single UAS, it is not possible to build a 3D model at a single instance of time since images are captured sequentially. As such, the 3D plume structure produced is only a time-averaged visualisation and some uncertainties at the margin of the 3D plume model are expected, due to turbulent movement of the plume during the period of imaging.
The estimation of 3D properties from sets of 2D images has gained more attention in recent years due to increasing usability of computer vision techniques. These include image segmentation and photogrammetric analysis. Stereoscopic and multiview photogrammetric analysis, which assess the movement of the projection plane when imaging from multiple angles, has gained attention for acquisition of 3D spatial properties of plumes. Examples include, using a combination SEVIRI and Meteosat-7 MVIRI data , and using MODIS and SEVIRI images for the Etna 2013 eruption . A combination of imagery from satellites in polar and geostationary orbits was used to determine cloud top height of the Eyjafjallajökull 2010 eruption (Zakšek et al., 2013).
The use of Structure from Motion (SfM) photogrammetry (Westoby et al., 2012) has also become increasingly popular due to the availability of commercial software packages with flexible implementations. SfM has been used extensively to analyse volcanic topography by generating Digital Elevation Models (DEMs) of volcanoes (James and Robson, 2012) (James and Varley, 2012) (Nakano et al., 2014) (Gomez, 2014) (Gomez et al., 2015) (Darmawan et al., 2017). Furthermore, SfM has been used for the proof-of-concept analyses of a forthcoming satellite mission, TOM (Telematics Earth Observation Mission), which is expected to retrieve accurate cloud top height measurements by acquisition of visible imagery from three different nanosatellites (Zakšek et al., 2018). A UASbased SfM method was utilized for reconstruction of a volcanic plume to retrieve its 3D properties (Gomez and Kennedy, 2018). In that study, the process was applied to a relatively stationary plume with low wind shear, using a single camera that was rapidly manoeuvred to minimize the evolution of the scene. The data in this work were acquired over a relatively long time period, yet it was still possible to estimate the 3D shape by assuming the plume occupied a similar space over the capture period. The method required interpolation of topographic data underneath the plume to give an indication of plume height as well as manual mesh manipulation. Although SfM is a powerful tool that can handle large number of views and missing data, the method is not always suitable for volcanic plumes which are often devoid of significant textures and have time-evolving features. Thus, since many plume datasets (including ours) are texture-poor, an alternative photogrammetric method is preferable. A space carving method, for building a 3D visualisation of volcanic ash plumes at Volcán de Fuego in Guatemala, was developed by the authors in previous work . In this previous study, simultaneous IR imagery captured from multiple ground-based cameras was used to estimate the plume's 3D convex-hull outline at any moment in time. From this convex-hull, quantitative measurements such as the speed and direction of plume transport could be deduced. This research investigates the application of this method to images captured from a UAS.

Regional setting
The images used for this study were acquired at Pacaya volcano, Guatemala (Fig. 1a). The modern edifice of Pacaya comprises a large summit crater that hosts the active MacKenney cone at its centre (Schaefer, 2012). After a long repose period between 1860 and 1960, intra-crateric Strombolian and effusive activity resumed in 1961 and has continued almost continuously since, interrupted by sporadic more intense explosive eruptions (Schaefer, 2012) (Rose et al., 2013). During image acquisition in February 2017, volcanic activity at Pacaya was characterised by persistent passive degassing and occasional weak explosions from the MacKenney cone. Extrusive activity from the base of the cone also fed several lava active flows that had begun to over spill the crater and descend the NW flank. As is common at arc volcanoes, the gas plume at Pacaya is composed predominately of water vapour (with minor transparent sulphur dioxide and carbon dioxide (Battaglia et al., 2018) and in condensed form this determines the visible margins of the plume ( Fig. 1(c)).

Space carving method
This section addresses the problem of creating a 3D model of a plume emitted from a volcano using sequential imagery acquired by a UAS-mounted camera. A process of space carving is implemented whereby a discretised search region is systematically classified by assessing the plume visibility in multiple segmented 2D images. This builds upon previous research (Wood et al., 2019) (with two important  The method for reconstructing the volcanic plume comprised pre-processing stages to transform camera pose (location and orientation) data and segment 2D images, followed by a two-stage 3D reconstruction. The first stage aims to locate the plume in the sky above the volcano, whilst the second stage performs detailed analysis to reconstruct the 3D model. differences: (1) each voxel in the scene is described by a probability value that classifies it as either a plume voxel or free-space, and (2) images were captured sequentially from a UAS-mounted camera. The workflow is shown in Fig. 2 with further explanation of the processing stages described below.
Space carving, also known as 'Shape from Silhouette' (Kutulakos and Seitz, 2000) is a generalised method that can be applied to reconstruct the 3D shape of an object using images captured from multiple known viewpoints. The method assumes that the boundary of an object in an image can be distinguished from the background to produce a silhouette image. It also requires that the camera locations and orientations are known a-priori. 2D segmentation and SfM can be used to achieve these requirements. The result of space carving is a visual-hull -the tightest bounding geometry of the shape of the 3D object and it is an outer-bound approximation of the object's shape. The principle is to create a 3D representation of an object, the visual hull, by sequentially projecting 3D space points into the silhouette images, then checking against a criterion.

Data acquisition
The image dataset was recorded on 21 February 2017 starting at 16:46 UTC for a period of 8-min and 10 s, with an image capture interval of 10 s. A total of 50 sequential images were captured as the UAS followed two successive circular flight paths centred summit crater at two different radii (42 m and 74 m) (see Fig. 1(d)). Although the volcano was continuously degassing during observations, the plume speed and transport direction remained stable and therefore the region occupied by the plume remained similar. A sample of the images acquired can be seen in Fig. 3.
Images were acquired using a DJI Phantom 3 Pro UAS featuring a built-in 12-megapixel DJI FC300X camera. This camera has a rolling shutter (Vautherin et al., 2016) which might introduce image distortion therefore reducing the overall accuracy. The effect is typically more detrimental when either items in the scene or the camera move with timescales similar to the shutter speed. In future, ideally a global shutter camera would be used. All images were geo-tagged with WGS84 coordinates at the time of capture using the UAS GPS.
A SfM software package (Agisoft Photoscan v1.4) was used to estimate camera orientations and optical properties. This was possible because the images not only captured the plume, but also the distinct terrain of the summit region. An additional benefit of using SfM was a further refinement of the camera positions therefore reducing potential inaccuracies introduced from GPS errors. The software outputs a file of xml format that includes the refined 3D camera poses and the camera's physical and optical properties.
The images were categorized into: (1) images showing the volcanic plume with no meteorological clouds present, Fig. 3(b), and.
(2) images showing the volcanic plume with meteorological clouds. The latter were further divided into two sub-categories: (2.1)images where the plume and background cloud appear largely connected (large overlap), Fig. 3(a) (2.2)images where the plume and background cloud are not attached to each other or only intersect slightly (small overlap), Fig. 3(c).
This categorization was used during the sensitivity study to describe the similarities and differences between four segmentation techniques (Section 3.1). Additionally, the first category served another purpose concerned with the identification of plume region prior to 3D reconstruction, as discussed in Section 3.2.

Coordinate transformations
Images were geo-tagged using latitude, longitude, and altitude information referenced to the WGS84 standard, however the space carving algorithm was implemented in a cartesian system, therefore the geometric transformations are used to convert between the two (Cai et al., 2011. The datum was the volcano vent, located at Latitude: 14.382°, Longitude: −90.601°, Altitude 2591 m. Once all image coordinates have been expressed in a single Cartesian system, known henceforth as the "East, North, Up" (ENU) system, sample points defined in this 3D 'world' space (X, Y, Z) can be projected into the camera coordinate system (x, y, z) using the camera GPS position and SfM orientation. This common 'world' space is a concept useful for later in the process. After this, points are mapped into the 2D image plane (x * , y * ) , and then quantised into pixel coordinates (u, v). This is the forward camera projection, and by assuming a pinhole camera model (Hartley and Zisserman, 2003), the sequence of transformations for the complete projection can be described using a mathematical model consisting of a set of matrix equations . Note that the pinhole model does not account for common distorting effects such as lens barrelling (Pears et al., 2012) (Luhmann et al., 2013). A representation of the various coordinate systems can be seen in Fig. 4.

2D image segmentation
Image segmentation is a process used to partition an image into multiple distinct regions or categories (Pal and Pal, 1993). Here it is used to identify the image pixels which contain plume. A good candidate image for segmentation contains regions of similar coloured pixels with sharp changes between regions. Several methods were trialled including: mean-gray level (Raju and Neelima, 2012), Otsu's method (Vala and Baxi, 2013), K-means clustering (El-Sayed, 2012), and semi-automatic. Mean-gray and Otsu's methods are automatic techniques based on image thresholding -that is, classifying pixels into two categories based upon their intensity value. Pixels with intensity values below the threshold are set to 0 (backgroundnon-plume) and those exceeding the threshold are set to 1 (foregroundplume). The result is a binary image; however, the contiguity of the segmented regions cannot be guaranteed. K-means clustering is also automatic and efficient in most cases and overcomes problems with non-contiguous foreground segments by introducing a pixel distance criterion.
The semi-automatic approach required user input to create binary images, hence, it enabled interactive refinement of the segmentation which could be advantageous, particularly when the background was not easily distinguishable from foreground and it consisted of many colours and high texture. MATLAB's (version R2017b) image processing software, 'Image Segmenter' (MATLAB, 2018) was utilized for semiautomatic segmentation. Binary images were created using tools such as, 'flood fill' and 'fill holes'. Using flood fill, the user only specified a starting point that represents a plume pixel and the method automatically selected patches with equal pixel intensity values. Additionally, a tolerance value could be specified based on Euclidian distance to control the extent of the selected area. The tool could be applied repeatedly until the user was satisfied the plume was completely segmented. The localised patchwork nature of this method eliminated the need to select a single global threshold value which can often result in mis-identified plume segments.

3D plume reconstruction
An algorithm based on space carving theory was developed for reconstruction of the 3D plume model using the UAS image set. Under ideal conditions, the object of interest would be timeinvariant, however due to the sequential nature of the UAS data, the plume shape evolves between each image. Because of this, the resulting 3D hull will be a time-averaged visualisation. Nevertheless, this assumption is deemed acceptable for this data set because the plume occupies roughly the same region during the entire data. For large and fast evolving eruptions, it might not be possible to sequentially image the plume using a single small UAS, under these conditions multiple much large UAS will become more appropriate to image the plume from higher altitudes and multiple directions simultaneously.
The method uses two stages of processing: The first is concerned with roughly identifying the plume region, and the second then focuses on a higher-resolution model within the sub-region identified in stage 1   5. A diagram representing the method of space carving. A region of 3D space (voxel) is projected into the 2D image (pixel) planes of cameras surrounding the plume. If pixel location is within the plume for all cameras, then the voxel is retained (green). If the pixel location is outside of the plume for one or more cameras, the voxel is discounted (red). The process is systematically repeated over a search volume, until only the voxels identified as plume remain . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (Fig. 2). During the first stage, the search region is defined by a coarse volumetric grid of 3D points (voxels) in the ENU coordinates, with the size of the voxel determining the resolution of the grid. It was assumed that the location of the plume was unknown, hence the initial search volume was set to cover a large region (see Fig. 6).
The centroid location of each voxel is projected into their pixel coordinates using camera projection. Then, each voxel is considered in turn and required to pass multiple criterions to qualify as a plume voxel. The first is the 'Voxel Visualisation Check'for a voxel, the algorithm examines whether the pixel coordinates lie inside image bounds. Images that cannot see the voxel are disregarded. Secondly, the 'Colour Consistency Check' is considered where the pixel coordinates of a voxel that meets the first criterion are checked if they reside inside the segmented silhouette of each image. Referring to Fig. 7, a white pixel is inside the silhouette and black is outside the silhouette.
For each voxel, a probability is calculated by dividing the number of silhouette images that recognise the voxel as a plume pixel over the total number of silhouette images that see the voxel i.e. for each voxel, it is the number of images that pass stage 2 (colour consistency check) divided by the number that pass stage 1 (visualisation check). This is a slight modification to the method described in  and was driven by the availability of many more viewpoints. Voxels below a probability threshold defined by the user (e.g. 75% for first iteration) are discarded, or 'carved' away, and only the plume voxels remain. The final stage is a 'Voxel Connectivity Check'if a voxel meets the threshold, but none of its neighbours are plume voxels, it is also discarded.
The first stage's output is a set of 3D points representing the detected plume region (Fig. 6). Once the plume region is located, the second stage follows the same process, but the new search volume is reduced to a smaller region which encapsulates the current detected plume plus a margin. This significantly reduces the computation workload. The limits for the new search region were defined by taking the maximum and minimum extent of the detected plume in the Cartesian axes. The second stage can then be processed using finer voxel resolution and more stringent probability threshold (e.g. 95%) resulting in a finer 3D model. The second stage voxel resolution was calculated using the relationship between the sensor pixel size, focal length and mean distance to object (plume), assuming a virtual image situated in the centre between the two circular flight paths.

Automatic and semi-automatic segmentation
A comparison between the segmentation results from different techniques is presented in Fig. 7. Four segmentation techniques were tested for two image categories; images where the plume region appears to merge with meteorological cloud (large overlap) and images where the plume region does not merge with meteorological cloud (small overlap) as per the description in Section 2.1. In all automated methods, regions with similar pixel intensity values to those of plume pixels, such as meteorological clouds or brightly exposed ground patches (brown regions in Fig. 7) are incorrectly classified as plume. Both meteorological cloud and plume were segmented as foreground when a large overlap (category 2.1) existed between the two regardless of the method used, whereas in cases of small overlap (category 2.2), the semiautomatic approach allowed a more precise selection of plume regions only as foreground disregarding the meteorological cloud or ground. Additionally, no ground patches were detected as part of the foreground, since the method incrementally segments regions of the image rather than applying a global threshold based on an intensity value. This shows the effectiveness of the semi-automatic approach (using human intuition) that enables successful segmentation of the plume. For repeatable monitoring, however, some greater level of autonomy could be introduced using more advanced image segmentation Fig. 6. Model setup showing the camera poses, the initial coarse search volume (25,972 voxels), the initially detected plume, and the refined secondary search volume (9,650,367 voxels). The coordinate system is ENU. algorithms. In the future a thermal camera might be beneficial to provide a temperature distinction between the plume and other image regions.

Plume identification
During the first iteration of the algorithm, the plume was successfully located using a coarse grid of voxels (20,000 to 60,000 voxels). However, this was achieved by eliminating images that incorporate meteorological clouds in the background (previously defined as category 2). In this case, 24 images (almost half the data set) contain meteorological clouds which have similar colour intensity to the plume. Hence, they pass both the colour consistency check and the 75% probability threshold, so the algorithm would have falsely considered meteorological clouds' voxels as plume voxels. With the meteorological clouds spanning a large region, the algorithm would locate the plume inefficiently narrowing only a small part of the search volume. The process for locating the plume region only required 10-20 s due to the low number of voxels and the reduced image set.
The first stage adaptive search led to a significant reduction of N70% in computational time for the second stage due to the smaller search volume. The identified search volume as seen in Fig. 6 has a volume of 2.49 × 10 5 m 3 . With a mean distance of 56 m, a focal length of 3.61 mm and a pixel size of. 1.5 μm, this corresponds to a resolution of 2.3 cm/pixel and 1.98 × 10 9 voxels. For comparison, performing detailed space carving on the original search grid (Fig. 6) of 1.87 × 10 6 m 3 volume directly would have increased the number of voxels to approximately 1.48 × 10 11 for the same spatial resolution. This caused the computational requirements to quickly exceed that available to most desktop users.

Final 3D plume model
The second stage of processing further refined the 3D model of the plume using a finer voxel search mesh. The algorithm was successful at determining the 3D structure of the emitted plume resulting in what can be effectively considered a point cloud of data (see Fig. 8); this can then be analysed further to retrieve its top height, volume, and drift direction (in this case southerly). It should be noted that these measurements are only applicable to the parts of the plume captured within the data; the whole plume is not captured in the source images. The top height of the eruption column reached 53± 4 m within the measured volume and is calculated by searching for the voxel with the greatest altitude above the summit. The plume volume was estimated to be 14,332± 560 m 3 , by multiplying the quantity of voxels with a voxel volume, that is -the voxel resolution to the cubic power. The plume shape is defined by the probability threshold and for 95% threshold value, the calculated plume distribution is presented in Fig. 8. Plume voxels with 100% probability were mainly concentrated in centre of the volume and further away from the vent. Since the plume starts narrow then disperses due to wind shear, a larger concentration of plume voxels exist as the plume moves away from the vent. Thus, there was a higher chance for the visual hulls of images to intersect and share common voxels further away from the vent leading to higher probabilities for those voxels to exist as plume voxels.
Using a Space Carving framework that assigns a probability for each voxel provided a more flexible approach than the original Space Carving framework (Kutulakos and Seitz, 2000). The original method did not incorporate a probability criterion, but rather required a simpler method whereby a voxel was classified as a plume if and only if all cameras that see the voxel agree that it is inside their silhouettes. As an example, if an arbitrary voxel (after projection to 2D) is perceived inside image bounds of 48 out of the 50 images, then using the probabilistic criterion, Fig. 7. Represents the Masked Images (Binary + RGB) using different segmentation techniques for two image categories of large and small overlap between moisture cloud and volcanic plume. Fig. 8. The plume voxels after the algorithm has been applied. The plume has been segmented along its midline (the red plane) to reveal the internal probability distribution. In general, more central voxels have higher probability. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the voxel is considered a plume voxel if it represents a foreground pixel in 46 binary images. For this example, the calculated voxel probability is 95.8% which is above the selected threshold of 95%. For the same example, the original framework would consider the voxel as a non-plume voxel because the projected voxel was outside just two of the silhouette images. The probability-based approach becomes more applicable with greater numbers of images, since with b10 it will become difficult to form a statistical significance. Setting the probability threshold to 100% allowed the retrieval of the results from the original framework as it implied that plume voxels were the ones where all corresponding visible cameras acknowledge them to be inside the silhouettes. The strict 100% probability threshold was less appealing because the plume did not have a defined shape due to its dynamic 'billowing' motion between image frames. This resulted in a less accurate representation of the plume.
In this study, it was not possible to build a 3D model at a single instance of time since images were captured sequentially. As such, the 3D plume structure is only a time-averaged visualisation and significant uncertainties at the margin of the 3D plume model would be expected, due to turbulent movement of the plume during the period of imaging.

Sensitivity to number of images
Studying the effects of number of images can be useful, particularly since it was recently shown that an increase in the image dataset can have detrimental consequences on photogrammetric analyses due to non-negligible evolution of the scene (Zakšek et al., 2018). The difference is that they were determining repeated surface models to track plume evolution; the work here incorporates equivalent variability through voxel probability.
To investigate the effect of the number of images used, two case studies investigating the consistency between the use of a subset of 35 and all 50 images is presented. The subset of 35 images were selected because they cover the first circular flight path ensuring that the set covers as many viewing angles as possible to minimize concavities and voids in the model. Fig. 9 shows the voxel clouds of Pacaya's plume using (a) 50 images and (b) 35 images respectively. The plume top height remained effectively the same at 53 m for both cases. In terms of plume volume, initially it was thought that the reduction in number of images would enlarge the plume shape, since the lower number of images implied a reduction in the number of intersections between the viewing cones generated by each image and hence, the visual hull of the object would become larger. However, a reduction of 254 m 3 in the plume volume when moving to 35 images was observed. This gave a similar indication to that noted by (Zakšek et al., 2018) where plume motion during image acquisition affected the optimum number of images used. Hence, a larger dataset may not always indicate better results. However, reducing the number of images further, caused an increasing estimation of the plume volume due to limited viewing angles.

Limitations and sources of error
The camera pose information will have a significant effect on the accuracy of the plume reconstruction. Initial camera positions were obtained using the GPS built into the UAS itself. Typically, these measurements might include multi-meter errors, however when compared to the size of the scene, in this case 100 s of meters, the error becomes less significant. Errors in the orientation will have a larger effect on the result. The orientations were estimated using a Structure from Motion (SfM) optimisation process and are potentially more accurate than any manual alignment. The SfM process could also be used to further refine the position estimates. It should also be noted that the 'plume top height' discussed in Section 3.5 represents only the upper bound at a certain distance from the ventit is not a height of equilibrium. The values given are defined by the image set at the particular time period and are given as an example of what might be possible with more expansive measurements.
In this work, each voxel was assigned a probability value that classifies it as either a plume voxel or a free-space voxel. However, the probability was driven by a consistency check, which decided how many images agreed that the pixel values corresponding to a voxel's projection lay inside the silhouettes (foreground). Hence, the accuracy of the result is significantly affected by this segmentation step and deteriorates for imprecise image segmentation. The quality of plume reconstruction is also impacted by the number of images used as discussed earlier. There is a benefit to be gained from more images at different viewpoints, however if these images are spread over a longer capture interval, the benefit might be outweighed by errors introduced due to the evolving scene. Hence, finding the optimum number of images to produce a close estimate of the true volume of the plume is still needed. Voxel resolution also has a large impact on the computational cost of the model given a fixed search volume. Errors will also arise from quantisation of the 3D space, but its significance will reduce as voxels get smaller. The use of additional optimisation methods, such as Octree data structure, could lead to an improved computational performance and these have previously been implemented in space carving algorithms (Montenegro et al., 2004) (Sainz et al., 2002). The adopted algorithm lacks information regarding plume motion, for example drift velocity or emission rates. Not accounting for plume motion in the analysis introduced an error that degraded the photogrammetric results which was unquantifiable at this stage. It was assumed that the plume broadly occupied the same region, hence any newly emitted plume would effectively replace the position of previous emissions. This is not completely true, since the effects of turbulence caused changes in the shape between images. In other volcanic scenarios, where the plume might be less persistent or visible, larger degradations are expected unless motion capture algorithms are employed. Optical flow has been used in previous work to measure factors such as turbulence and emission rates (Thomas and Prata, 2018) (Gliβ et al., 2018) (Peters and Oppenheimer, 2018), and extensions to that work might enable similar analysis in 3D. To completely remove the errors due to plume evolution, the images should be captured simultaneously. This could be achieved using multiple synchronised UAS, however the field implementation of such a system would be challenging at this time.

Conclusions
A 3D volume estimation algorithm has been developed based on an image dataset capturing a volcanic gas plume emitted from Pacaya volcano, Guatemala. The methodology has been presented, including full details of the coordinate transformations, camera projection, image segmentation, and 3D reconstruction steps. At its core, the algorithm developed used a previously implemented space carving framework , but with several important modifications; specifically, implementing a probability-based plume identification criterion, and a two-stage localisation refinement. The method relies on image segmentation and multi-view imagery for reconstruction, with the computational cost based on the specified voxel resolution, number of images, and plume search region.
Applying the algorithm to images of the volcanic plume from Volcán Pacaya enabled determination of the plume shape and physical properties, including top height, drift direction, and volume. Of the observed sub-section of plume extending approximately 100 m from the vent, the top height reached 53 ± 4 m, with a southerly drift, and a voxelbased volume of 14,332± 560 m 3 . These properties are derived directly in 3D rather than relying on approximations of the plume structure that are often applied to 2D analyses. The results presented confirm that the 3D topology of volcanic plumes can be measured using single UASbased photogrammetry, assuming the plume motion is persistent and broadly occupies the same region throughout acquisition. Uncertainties related to time-averaging of the plume visualisation are minimized further by initializing the algorithm with a large dataset of images at multiple viewing angles.
The adopted method has direct applications in volcano monitoring since it presents a novel approach to constrain volcanic source parameters. The deployment of multiple UAS may allow the retrieval of timevarying 3D spatial properties. With this technique advancement, the plume 3D convex hull could be calculated for each time-step, allowing the plume evolution and dynamic behaviour to be studied. UAS-based photogrammetry, and the potential future developments within this field, may contribute crucial data with which to validate plume dispersal models.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.