Monitoring and Computation of the Volumes of Stockpiles of Bulk Material by Means of UAV Photogrammetric Surveying

Abstract

The monitoring and metric assessment of piles of natural or man-made materials plays a fundamental role in the production and management processes of multiple activities. Over time, the monitoring techniques have undergone an evolution linked to the progress of measure and data processing techniques; starting from classic topography to global navigation satellite system (GNSS) technologies up to the current survey systems like laser scanner and close-range photogrammetry. Last-generation 3D data management software allow for the processing of increasingly truer high-resolution 3D models. This study shows the results of a test for the monitoring and computing of stockpile volumes of material coming from the differentiated waste collection inserted in the recycling chain, performed by means of an unmanned aerial vehicle (UAV) photogrammetric survey and the generation of 3D models starting from point clouds. The test was carried out with two UAV flight sessions, with vertical and oblique camera configurations, and using a terrestrial laser scanner for measuring the ground control points and as ground truth for testing the two survey configurations. The computations of the volumes were carried out using two software and comparisons were made both with reference to the different survey configurations and to the computation software.

1. Introduction

A fundamental process in many activities is the management, handling and storage of bulk goods of various kinds. This can be seen in many business fields in which bulk materials constitute a cost, a resource, or is a phase of a process, e.g., in civil engineering [1], mining and extraction [2], industry [3], agriculture, food [4] and in the management of waste landfill and recycling sites [5].

In the past, the continuous monitoring and computing of the volumes was conducted with classical topography tools using Total Stations or global navigation satellite system real-time kinematic (GNSS-RTK) systems. These techniques are still used today but pose problems related to the accessibility and security of the sites, to the survey timings and relative costs. Moreover, in this way it is possible to measure few characteristic points only, not enough for adequately describing the most complex shapes.

Section 2 examines the survey of storage sites of bulk materials and their computation methods. The laser scanner and photogrammetric techniques (terrestrial or UAV), now widely used in many fields [6,7,8], instead make it possible to carry out rapid and dense samplings of complex surfaces without the need to access directly to the sites, cutting operational costs and ensuring increasingly accurate models and greater safety. Consequently, these techniques are now extensively replacing the traditional ones, also in the context of continuous monitoring [9] and computation of the stockpile volumes of accumulated or excavated materials and cut and fill backfilling [10].

Section 3 presents the current case study, related on a feasibility test of the use of UAV photogrammetry for the frequent monitoring of the volumes of waste deriving from differentiated collection accumulated in a storage and pre-processing site for recycling (Figure 1). Section 4 shows the results obtained with different camera networks (vertical, oblique and vertical + oblique acquisitions) and volume computation methodologies. In Section 5, the results are compared with each other and with terrestrial laser scanner acquisitions and assessed in terms of timing, reliability and performance for determining the most effective workflow.

2. Survey and Computation Methods

Irrespective of the techniques used, the main steps that lead to determining the volume of piles of bulk materials include:

• the surveying of the surfaces through acquisition of points and characteristic lines;
• the modelling of a solid that approximates its real shape and size;
• the geometrization and discretization of the solid into elementary volumes.

The methods via which each of the three phases are carried out highly depend on the purpose of the computation and on the constraints placed by the form, size, position and possible ways of storing the volumes themselves. As a result, it is important to carry out a careful preliminary analysis and a survey of the locations in the aim of evaluating the possible options and appropriate methodologies.

The accuracy of the volume computation will be influenced by:

• the number and accuracy of the acquired points;
• the method of modelling the surfaces;
• the knowledge and identification of the underlying ground;
• the shape and complexity of the stockpile surfaces.

2.1. Survey Methods

In the past, the survey of accumulated materials was carried out using traditional topographic methods, and more recently via use of Total Stations and GNSS-RTK technology. They allow accurate and fast measures of few significant points selected by the operator, so are cost-effective for measuring simple piles or small sites. In more complex cases, a great amount of points has to be measured to reconstruct the surface and measure the volume of a stockpile, requiring a longer time and higher costs. Moreover, sometimes it is difficult or unsafe to acquire all the needed points because the operator has only a partial view of the site from the ground, and it is not possible to go up on the piles because they are steep and unstable and survey operations interfere with other production activities.

Laser scanner and close-range photogrammetry, today widely used in many survey application (land, infrastructure, architecture, archaeology and so on), remotely and rapidly measure high density indiscriminate points. The increasing hardware computing power and high-performing software allow the processing of large volumes of data.

Scanning systems based on LIDAR (Light Detection and Ranging) technology acquire per unit of time a large number of 3D points on objects of any size and complexity, with very high resolution and great accuracy [11,12]; each scan directly creates a point model (also called a point cloud) that must be usually aligned with other scans for generating a complete 3D model (Figure 2 and Figure 3).

Additionally, close-range photogrammetry is now widely used for rapid, accurate and low-cost sampling of complex surfaces, irrespective of the extent of the area acquired [13]. Starting from digital imagery (even from non-metric cameras) and thanks to Structure-from-Motion technique [14,15], the internal and external orientations of the acquired images are reconstructed via the recognition of thousands of homologous points in the overlap areas of the images. Like in laser scanning, dense matching algorithms produce point clouds with highly automated and geometrically robust procedures.

These models are scaled and referred via GCPs (Ground Control Points) acquired in a pre-established reference system. Several known coordinate points, not used for the computation of the orientation and the scaling, allow for assessment of the accuracy [16].

Both laser scanner and photogrammetry can be used from static ground stations (sometimes permanent for continuous monitoring), or in mobile or airborne configurations (on conventional or remote piloted aircraft), depending on the purposes of the survey and the operational, topographical and morphological conditions of the sites.

For the monitoring and volume computation of piles of bulk materials, the two techniques have been used both from the ground and on UAVs. Respect to the above-mentioned Total Station and GNSS survey, they have been preferred because it is possible to get a denser sampling and therefore a more reliable volume computation, even of the most complex and wide area. As laser scanners and photogrammetry acquire only in the line of sight, in many sites it is not possible to obtain a complete volume reconstruction from the ground only. For this reason, UAVs are always more frequently used for acquiring stockpile storage areas from above. As most cameras are lighter than laser scanners, UAV photogrammetry is more affordable than laser scanning because the lower payload makes it possible to use low-cost UAVs (Figure 4).

Numerous experiences document the best results obtainable in terms of accuracy of the volumes computed with close-range photogrammetric surveys compared to total station or GNSS [17,18,19], as well as evident advantages also in terms of economic and time costs [20,21] and the achieving of higher levels of reliability and repeatability [22] than with traditional techniques.

In any case, choice of the most appropriate technique depends on the purpose of the survey, the dimension and the morphological and topographical conditions of the site, accessibility, visibility, safety and a careful cost-benefit analysis. In many cases, surveys integrate different methodologies and GNSS and Total Station are needed to measure control networks, Ground Control Points or Check Points for assessing, scaling and georeferencing the acquired data.

2.2. Volume Computing Methods

Three different methods for computing volumes are usually considered [23]:

• by cross section;
• by horizontal section;
• by prisms between two surfaces;

where the choice among them depends on specific solid.

The cross-section method, typically used for fairly regular solids, computes the volume between two cross-sections of an established area. The computation method may differ depending on the characteristics of the solid, according to more or less complex formulas of the following type:

$$Vi=\frac{d_1}{6}\left(P_1+4{\left(\frac{\sqrt{P_1}+\sqrt{P_2}}{2}\right)}^2+P_2\right)$$

(1)

or, in simpler cases:

$$Vi=d_1\frac{P_1+P_2}{2}$$

(2)

with reference to Figure 5.

The horizontal section method uses the classical contour lines placed at an appropriate equidistance; the volume is divided into more elementary solids by horizontal planes corresponding to the contour lines (Figure 6). The volume is computed with formulas similar to Equations (1) or (2), depending on the morphology of the solid.

The method via prisms between two surfaces consists of the definition of two surfaces, determined by means of interpolation of 3D spatial data, one of which has the function of real surface (upper surface) and the other of a reference surface (lower surface); the solid delimited by these surfaces is discretized into elementary prisms with a triangular base (Figure 7) or trapezoidal prisms from which volume is computed using the following formula (case of a triangular prism):

$$Vi=P\frac{{{\displaystyle \sum }}_{i=1}^{3}Hi}{3}$$

(3)

For objects with a predominantly linear development (e.g., road prisms), the cross-sectional method is commonly used. For trapezoid-shaped objects, both the trapezoidal and cross-sectional methods can be used.

The method via prisms between two surfaces is the one that will be discussed below and which is the basis of the computation principles of the software used and described in this paper.

2.2.1. The Influence of the Surveying Method

The surveying method affects both the amount of acquired points and their positional accuracy; the reliability of the resulting surface model to the real conditions depends on the combination of these two factors.

With reference to the number of points acquired, different experiences clearly show that the best results are achieved with laser scanner [24,25], or close-range ground [20,26] or UAV-based [27] photogrammetric survey techniques, which allow for acquiring a considerable density of points and generating models that are closer to reality. This means that the volume computing procedures give more accurate results.

2.2.2. Elevation Models

Among the approaches available for digital terrain modelling (DTM), the two most frequently used for model computation with the prism between two surfaces method are:

• vector-based data modelling (Triangulated irregular network - TIN);
• raster-based data modelling (Raster Digital Elevation Modelling - DEM).

TIN modelling is generally more flexible as it can better describe break-lines and abrupt slope variations, thus allowing for a more precise result. It is created starting from the original points, except in the case of subsequent simplifications and interpolations, and each point is the vertex of a triangle.

In the raster-based regular grid form (Raster DEM, see Figure 8), the model has no directly measured points, due in any case to being interpolated, and it is a less flexible model in representing sudden changes of slope. The dimension of the grid influences the accuracy of the model itself; the smaller the cell, within certain limits, the more representative the model will be. The dimension must nevertheless be related to the aim of the survey, the complexity of the surfaces to be described and, last but not least, the amount of data: the smaller the grid dimension, the more data will need to be processed and the slower the processes will be.

2.2.3. The Definition of the Reference Base

As described, the computing method for prisms between two surfaces provides for the definition of an upper surface and a lower surface; if the upper surface is in any case easily determined with the above-mentioned methodologies, a critical point represents the definition of the lower surface, i.e., the reference base. In fact, this is not directly detectable and, unless it is known a priori, it must be determined by interpolating significant points of the uncovered ground surrounding the volume to be computed. While on one hand it gives good results for relatively flat base surfaces, this operation can be really challenging for irregular surfaces with sudden changes in altitude or slope; this criticality is further accentuated by the fact that an elevation uncertainty at the base of the stockpile produces much greater errors in the volume estimation of the same uncertainty concentrated at the top [28].

The determining of the base surface is therefore an aspect that should not be underestimated and which must be taken into account in the overall survey project and computation flow. For regular surfaces very close to a flat or to a rolling plane, the currently most used solutions provide for performing a high-precision survey of the most significant points of the surrounding uncovered surfaces and through these creating interpolated ground models (TIN, Raster DEM) [22].

In the case of a general survey performed with laser-scanner or photogrammetric techniques, the points necessary for determining the reference base can be measured by means of a specific GNSS-RTK survey, otherwise they can be selected in the area surrounding the stored material in the same point clouds acquired with the techniques used for surveying the overall area.

2.2.4. The Influence of the Shapes and Size

In addition to the factors described above, the accuracy of the final computation of the volume depends on the shape and size characteristics of the solid, therefore it is not possible to determine in absolute terms but only in relative terms. The accuracy in the volume computation depends on the ratio between surface and volume [29,30]: the smaller the ratio, the more accurate the volume (Figure 9).

3. Application Experience in a Stockpile Area of Recyclable Waste

3.1. Foreword

This study has been conducted as a preliminary feasibility test of a permanent monitoring project, to be repeated monthly, for determining the variation in quantity of the material stored in a recyclable waste collecting, processing and storage plant.

The site has a total area of about 15,000 square meters, with materials stored both indoors and outdoors (Figure 10). The survey involved stockpiles of materials located over an uncovered area of about 10,000 square meters on a reinforced concrete floor with a regular and practically flat surface, except for slopes for draining the wastewater. The stockpiles measured consist of glass waste already pre-processed and reduced to fragments. The material is opaque and for the data acquisition purposes is comparable to sand or gravel with a granulometry ranging from less than a millimetre to a few centimetres. Part of the stocking area (visible in the upper left corner of Figure 10) is covered with sheets of plastic material that reflect the light. The volumes to be computed are partly delimited by permanent structures (fixed retaining walls) or provisional barriers (delimitation with cubes in concrete or other material), and partly piled against the boundary wall of the property. From the site configuration described, it is evident that the volumes are not accessible 360° for a complete survey from the ground, however they are clearly visible from above. In addition, the plant is continually affected by the circulation of workers and vehicles for the mechanical loading/unloading and handling, so a survey carried out at ground level would be both obstructive and hazardous.

All the above-mentioned site characteristics led to select UAV photogrammetry for a feasibility test. In the most recent studies on the UAV photogrammetry applied to the topographic survey of quarries [31] or geomorphic features, in addition to the traditional configurations with vertical camera have sometimes been used also inclined cameras. In these cases, it improved the representation of slope geometries and therefore it was decided to test if there were significant advantages in the stockpile’s measurement also.

Every photogrammetric project requires to be scaled and validated inserting a certain number of GCPs (Ground Control Points) and CP (Check points), to be measured with other technique. For a permanent monitoring activity, the best solution would be to set up a stabilized geodetic network constituted by durable targets uniformly distributed over the area, measured (and periodically checked) with total station or GNSS.

As this was a feasibility test, in view of a future monitoring implementation on a regular basis, it has been considered sufficient for first evaluations to use as GCPs the coordinates the centre of temporary targets acquired with a concurrent laser scanner survey. This has been also useful for conducting tests and comparisons between methodologies and for evaluating the accuracy achieved. The uncertainty of GCPs coordinates acquired from laser scans is higher than if they were measured by total station or GNSS, but all targets were acquired from a short distance and, above all, all the following calculations use the same coordinates and are consequently affected by the same systematic error. The inaccessibility of some areas of the site for a complete terrestrial laser-scanner survey made it impossible to place the targets uniformly over the entire area and therefore some areas of the photogrammetric model are not properly covered by GCPs; nevertheless, the GCPs configuration implemented was deemed adequate for the purposes of the test.

3.2. The Survey

The survey was performed on the same date both with terrestrial laser scanner and the UAV photogrammetry. The survey date was chosen on a summer Saturday in order to perform operations in conditions of minor interference with the working activities of the site, and the acquisition lasted around four hours.

Beforehand, nine temporary targets, A3 format, black and white, arranged in a horizontal position on the ground or on the boundary walls, were distributed over the area to be surveyed, functional both to the laser scanner survey for the following alignment of the scans, and the photogrammetric survey as Ground Control Points and CPs (Figure 10).

The laser-scanner survey was performed with a Zoller + Fröhlich 5010C scanner [32]; the number and position of the scans was determined according to the morphology of the terrain, the geometry of the area, and the limited visibility of some of the surfaces (Figure 11). Due to the lack of suitable vantage points, 18 scans were requested for the complete acquisition of the area, four of them with a sampling step of 6 mm at 10 m, and the others of 3 mm at 10 m, only acquiring the 3D coordinates and the reflectance value, as RGB value was not useful for the aim of the survey. With a view to aligning the scans mainly using natural points, the lower resolution scans were carried out in some areas non-relevant for the survey but to ensure the continuity of the survey.

The summary data of the survey are listed in

Table 1. Laser scanner survey data.

Data	Quantity
Number of stations	18
Nominal resolution	no. 14 SH (3 mm/10 m) no. 4 H (6 mm/10 m)
No. of operators	2
Survey duration	Approx. 3 h

.

For the photogrammetric survey, the DJI Phatom4Pro multirotor quadcopter UAV [33] was used and the flight plan was created on the UgCS software [34], setting the parameters that allowed for acquisition of images with a GSD (Ground Sampling Distance) of approximately 2 cm (Figure 12). Two series of images acquisition were programmed and carried out, one with the camera in a vertical position and the other with a 30° angle.

Both acquisition phases took approximately four hours, including the preparation, test, battery recharging and quality control of the acquired images. The preparation phase also entailed the identification and preparation of a ground control station that would not interfere with the ordinary activities of the storage site. The duration of the flights was about 15 min each.

Table 2. Parameters of the photogrammetric acquisition.

Data	Quantity
Flight altitude	Approx. 75 m from ground level
Forward overlap	85%
Side overlap	70%
Number of images	139 (vertical camera) 144 (30° oblique camera)
Duration	15 min. per flight
Survey duration	Approx. 4 h

contains the main parameters of the flights planning and execution.

3.3. Processing

3.3.1. Laser Scanner Survey

The laser scanner acquisitions were imported and processed with Leica Cyclone software [35]. The scans were aligned in a local reference system chosen at the origin of Scan No. 4 because it was approximatively in the centre of the area of interest, and both natural points and the nine targets were used for this activity (Figure 13); the final model is composed of 1,308,042,626 points. The scans registration “mean absolute error” is 1 mm (Figure 14).

From the point model the 3D coordinates of the nine targets were obtained for being used as known coordinate points (GCPs) to scale and orientate the photogrammetric model. Considering the registration error, the scanner characteristics and the scanner-target distance, GCPs have a sub-centrimetric accuracy.

Finally, a list of the coordinates of the targets was produced in the local reference system for being imported directly into the photogrammetric image processing software, in this way all 3D data, obtained both from laser scanning and photogrammetry, are in the same reference system and can be directly compared and integrated.

3.3.2. Photogrammetric Survey

The processing was carried out with Photoscan software [36]; all the phases described below were performed for each of the following three scenarios:

• Vertical camera images only;
• Oblique camera images only;
• Vertical + oblique camera images.

After being imported into the software, the images were examined to make a rough evaluation of their quality, eliminate any unsuitable ones, and carry out some masking interventions on possible disturbing elements. After image alignment, GCPs have been collimated on images imposing them the coordinates determined in the local reference system of the laser scanner survey, and by means of which it was possible to achieve the absolute external orientation and scaling of the model.

Of the nine collimated targets, only six were used as GCPs, while three of these were excluded from the calculation and used as CPs (Check Points). The statistical values of the external orientation operation reported a total RMSE value of 0.54 cm on the GCPs, and 1.60 cm on the CPs for the vertical flight; the target 05, taken as a CP, was in a hardy visible position, and the major errors were noted on this one.

The following

Table 3. GCP errors in the vertical-only acquisition (A).

GCP	X err. (cm)	Y err. (cm)	Z err. (cm)	Tot. err. (cm)	Projections	Image (pix)
01	−0.011	−0.253	−0.300	0.393	45	0.093
03	−0.244	−0.176	−0.790	0.845	52	0.155
04	−0.121	0.163	0.121	0.236	33	0.102
06	0.165	0.100	0.156	0.248	49	0.116
08	−0.103	0.023	0.021	0.107	45	0.100
09	0.321	0.146	0.812	0.885	48	0.142
Total	0.189	0.160	0.485	0.545		0.122

and

Table 4. Check point errors in the vertical-only acquisition (A).

CP	X err. (cm)	Y err. (cm)	Z err. (cm)	Tot. err. (cm)	Projections	Image (pix)
02	0.110	0.794	0.657	1.036	36	0.106
05	0.843	−0.726	−2.278	2.535	48	0.140
07	0.213	0.341	0.172	0.437	51	0.128
Total	0.506	0.652	1.372	1.601		0.127

show the error parameters relating to the GCPs and CPs for the vertical flight (A), the errors for the oblique and vertical + oblique flights (B and C, not reported) have the same order of magnitude.

For each of the three scenarios, the obtained dense clouds have a density of approximately 1000 points/square meter (Figure 15).

Each dense cloud was edited to eliminate points belonging to not relevant objects that could have influenced the volume computation, for example the branches of trees growing on the border (Figure 16).

Finally, starting from the dense clouds for each of the three scenarios, the DEM of the entire surveyed area was generated, necessary for the subsequent computation of the volumes.

The DEM of the reference base was generated by exploiting the same dense clouds created previously; starting from these, a manual classification of the terrain was carried out by selecting the areas interactively recognised as such and then creating an interpolated 3D model only using the points classified as “ground” (Figure 17).

At the time of the survey, the areas not covered by stockpiles were reduced to a minimum due to the effect of a significant contribution of typical material during the summer. The repeating of the process in a period of minimum load could have a broader terrain classification and therefore ensure greater fidelity with the real conditions. In any case, given the conditions of substantial flatness of the terrain, the solution, even in the conditions at the date of the survey, was satisfactory.

At the end of the processing phases, the two surfaces necessary for the computation (upper and lower surfaces) were generated in the form of a DEM grid with the same resolution.

3.4. Volume Measurement

The volumes were calculated for each of the three scenarios with two different software:

• Agisoft Photoscan;
• Esri ArcGIS [37].

Both packages require the two computed DEMs (in Photoscan the ground one is necessary for determining an average ground level) as well as for determining the volumes to be calculated. The volumes were delimited in Photoscan in such a way as to use the advantages of a 3D acquisition directly on the model generated and make sure that the vertices lie on the volumes to be calculated; the boundaries were subsequently exported as shapefiles to be reused in Photoscan itself as well as in ArcGIS.

All in all, nine areas were extracted, in addition to constituting a positive and complete test base, the subdivision was also necessary for differentiating the areas based on the type of material stocked because the final objective was that of computing the volumes differentiated depending on the characterization of the waste products.

3.4.1. Volume Measurement Procedure with Photoscan

In Photoscan, the “reference surface” to be used for the computation of the volumes, can only be a plane and the insertion of a more complex reference surface determined a priori is not envisaged. This does not constitute a limitation for situations, such as in this case study, where the base surface is effectively flat or extremely similar to a plane (Figure 18).

This flat surface can be defined in three ways:

• as an inclined flat surface and interpolated by the vertices of the delimitation of the stockpile (best fit plane);
• as a horizontal flat surface with a medium height determined by the heights of the vertices of the delimitation of the stockpile (mean level plane);
• as a horizontal flat surface at a reference height defined by the user (custom level plane).

The first two options are only well suited to stockpiles for which the base contour is totally acquirable and therefore not delimited by any containment artefacts. However, due to the fact that all the stockpiles on the site are in some way delimited by some sort of artefact or other, the third option was used (custom level plane).

Using the “custom level plane” calculation option, the value entered was the average value of the base DEM on the single delimitation; this value, which cannot be calculated directly in Photoscan, was computed for all the boundaries with external GIS software. The volumes of the nine areas of interest were interactively computed in this way (Figure 19).

3.4.2. Volume Measurement Procedure with ArcGIS

In ArcGIS, the measurement was made using the “cut fill” tool. The principle is the one already referred to above and used to determine the volume between two surfaces (DEM raster) through the division of the volume into elementary prisms, the base of which is defined by the resolution step of the DEM, and the height of which is equal to the difference between the altitude information of the two DEMs in the same cell (Figure 20).

In ArcGIS, the processing steps of Photoscan were therefore imported, namely:

• the boundaries, contained in a single shapefile with nine area geometry records;
• the DEM of all the surfaces;
• the DEM of the reference base.

To perform the computation on the nine areas of interest, a workflow was created using the “model builder” tool which, by means of iteration, performs the extraction (clip) of the surfaces for each area of interest and carries out the computation by generating a general summary table (Figure 21). This guarantees having an automated and repeatable procedure at the end of the processing in Photoscan that supply the input data, while limiting the interactive procedures to a maximum.

4. Results

Volume Assessment

The values calculated were compared both for assessing the differences in results between the two software types used for the same scenario, and for evaluating the differences in the values calculated between the three acquisition scenarios.

As regards the comparison between the types of software used, we took the calculation with ArcGIS as a reference, and therefore, the positive differences indicate a greater volume computed with Photoscan. A substantial equivalence of results was obtained between the two software, the maximum absolute differences were found for the volume “6” between +2.40 m³, and +3.05 m³ (equal to a percentage difference between 0.04% and 0.05%). The maximum percentage differences are all below 0.2% (0.17% for volume “4” in the vertical + oblique scenario); these differences, which are not statistically significant, lead to an evaluation of the results obtained with the two software as effectively the same. This also gave rise to the evaluation of the simplification performed in Photoscan (by imposing a flat base surface with a reference level equal to the average level of the surface model of the base) as not producing any appreciable effects in the case of ground surfaces that were actually almost flat as the one in question.

Table 5. ArcGIS vs. Photoscan volume computation in the three flight configurations.

	A—Vertical Camera only				B—Oblique Camera only				C—Vertical + Oblique Camera
ID	ArcGIS m3	Photoscan m3	Diff. m3	%	ArcGIS m3	Photoscan m3	Diff. m3	%	ArcGIS m3	Photoscan m3	Diff. m3	%
0	394.87	394.96	0.09	0.02	397.12	397.07	−0.05	−0.01	398.40	398.51	0.11	0.03
1	460.56	460.46	−0.10	−0.02	463.12	463.10	−0.02	0.00	465.05	464.94	−0.11	−0.02
2	388.04	388.09	0.05	0.01	392.29	392.38	0.09	0.02	391.09	391.21	0.12	0.03
3	1247.64	1248.10	0.46	0.04	1253.00	1253.56	0.56	0.04	1256.03	1256.49	0.46	0.04
4	77.29	77.33	0.04	0.05	79.61	79.63	0.02	0.02	74.75	74.88	0.13	0.17
5	616.58	616.71	0.13	0.02	627.99	628.34	0.35	0.06	619.57	619.71	0.14	0.02
6	5512.90	5515.72	2.82	0.05	5571.03	5573.43	2.40	0.04	5559.63	5562.68	3.05	0.05
7	2049.24	2049.14	−0.10	0.00	2074.82	2074.80	−0.02	0.00	2060.83	2060.95	0.12	0.01
8	1519.61	1521.93	2.32	0.15	1532.28	1533.82	1.54	0.10	1535.67	1537.01	1.34	0.09

summarizes the volumes calculated with the absolute differences and percentages among the values obtained with the two software.

As far as the differences between the results obtained in the three scenarios are concerned, given the perfect equivalence between the used software, for the sake of brevity we will only refer to the values calculated with ArcGIS.

The values obtained with the vertical acquisition scenario were used as a reference; therefore, positive differences refer to higher volumes in the other scenarios.

The absolute differences in the module are found between 2.25 m³ (between oblique and vertical acquisition in volume “0”) and 58.33 m³ (between oblique and vertical acquisition in volume “6”), with a tendency, although not always respected, towards greater differences for the larger volumes. These differences, which may seem relevant, above all in the maximum values, are however much less significant when expressed in percentage terms.

Percentage differences vary in the module between 0.43% (between oblique and vertical acquisition in volume “3”) and 3.29% (between vertical and oblique acquisition and vertical acquisition in volume “4”), which is consistent or even better than what usually results from the literature for similar cases with regard to the accuracy of the method [17,38,39].

It should be noted that the highest module value is found in volume “4” which, compared to the others, has a greater ratio between surface and volume (about double) which is in line with what has been stated in Section 2.2.4 about the correlation between accuracy and surface/volume ratio.

Table 6. Comparison among the volumes computed in the three scenarios.

		A Vertical	B Oblique	C ver. + obl.	B-A		C-A
ID	Area m2	Volume m3	Volume m3	Volume m3	Volume m3	%	Volume m3	%
0	197	394.87	397.12	398.40	2.25	0.57	3.53	0.89
1	221	460.56	463.12	465.05	2.56	0.56	4.49	0.97
2	185	388.04	392.29	391.09	4.24	1.09	3.05	0.79
3	457	1247.64	1253.00	1256.03	5.36	0.43	8.39	0.67
4	86	77.29	79.61	74.75	2.32	3.00	−2.54	−3.29
5	267	616.58	627.99	619.57	11.40	1.85	2.99	0.48
6	1401	5512.90	5571.03	5559.63	58.13	1.05	46.73	0.85
7	461	2049.24	2074.82	2060.83	25.58	1.25	11.59	0.57
8	753	1519.61	1532.28	1535.67	12.67	0.83	16.05	1.06

contains a summary of areas of the piles, their volumes (as indicated in Table 5) and the differences found among the three scenarios, with the minimum values (in green) and maximum values (in red) highlighted for each calculation, while the values of the surface/volume ratio are illustrated in

Table 7. Surface/volume ratios of the nine stockpiles (in C scenario).

S/V Ratio
0	1	2	3	4	5	6	7	8
0.49	0.48	0.47	0.36	1.15	0.43	0.25	0.22	0.49

.

5. Discussion

5.1. Comparison Among the Surface Models

To more effectively analyse the differences, certainly attributable to the surface models, comparisons have been carried out in various ways on the point clouds and DEM produced.

Several vertical profiles have been created both on the point clouds and on the DEM and the comparisons carried out; where possible they have also been compared with the profiles coming from the laser-scanner point cloud considered as ground truth.

Initial qualitative considerations have been carried out comparing the point clouds generated in the A, B, and C scenarios with the laser scanner point cloud (Figure 22). In comparison with the model of laser-scanner points, both for the clouds of the vertical acquisition (A) and for those of the vertical + oblique acquisition (C), average differences were found (excluding outliers) between 2 cm and 3 cm (in absolute values) in practically the entire area with a slight preponderance for negative values; said values only increase on average to 5–7 cm for the oblique acquisition (B).

Additionally, in reference to the laser scanner point cloud, we instead found average differences of about −7 cm for the vertical acquisition (A) and around +7 cm for the vertical + oblique acquisition (C) localized at the top of volume “6” (which shows much greater differences in m³), with a consequent and decidedly elevated difference of about 14 cm between the two scenarios A and C (Figure 23).

This difference is confirmed by comparing the profiles of the DEM generated with Photoscan in the three scenarios.

In the remaining areas, the differences between the photogrammetric DEMs are maintained in the order of 2–3 cm. The profiles also show a strong noise component at the top of volume “6” in the DEM of the oblique only acquisition (B, see Figure 24).

From the general considerations, as a first result has been considered that the oblique only acquisition (B) gave worse results, so the further investigations described below were conducted on the vertical acquisition (A) and vertical + oblique acquisition (C) DEMs only.

In particular, in order to support the results of the comparison of the profiles with statistically rigorous data, comparisons were carried out with the CloudCompare software [40] directly among the point clouds generated with the photogrammetric shots (Figure 25).

For the comparison, only the areas involved in the volume computations were selected, and the areas affected by the stockpiles from nos. “0” to “7” (contiguous) were separated from the stockpile area “8” (isolated).

The comparison was made using the M3C2 plugin [41], which computes the signed distance directly from two-point clouds, by indicating the cloud of the vertical acquisition as a “reference” and the cloud of the vertical + oblique acquisition as “compared”. In the comparison of the area involved with the stockpiles from “0” to “7”, approximately 4,000,000 points were compared with an average value of the differences in Z of 3.9 cm and a standard deviation of 4.2 cm.

When comparing the area affected by stockpile “8”, approximately 1,000,000 points were compared with an average value of the differences in Z of −0.5 cm, and a standard deviation of 1.8 cm, which are decidedly better than the previous ones (Figure 26).

The spatial distribution of the waste products found for stockpiles from “0” to “7” confirms what was already found by characterizing the profiles, that is, there were larger and more significant waste products located on the top of stockpile “6” with values of between 10 cm and 20 cm.

With the exception of the “6” stockpile area, there is consequently a good compliance between the two clouds with statistically non-significant differences. The excessive differences located in this area repeatedly highlighted are plausibly due to a lack of GCPs nearby and therefore to a certain effect of model drift; this problem should not arise in the possible future application of the proposed method for the periodic monitoring of the area, in which case the positioning of the GCPs would be made more homogeneously with a uniform distribution over the whole area. In addition, as already mentioned, the area is partly covered by an optically non-cooperative tarpaulin, the reflections of which could have had an impact on the results obtained both with the photogrammetry and the laser scanner. From the outcome of the comparison it can already be affirmed that the two acquisition configurations A and C are to be considered qualitatively equivalent; in order to better support this assessment and above all, for the purpose of evaluating errors on the calculated volumes, the same criterion was used to compare the points model obtained by the laser scanner.

5.2. Estimation of the Errors in the Volume Computation

As can easily be intuited, and as demonstrated by several studies [38,39], the estimation of the uncertainty in the volume computation can be correlated with the evaluation of the uncertainty of the altimetric models generated by the surveying and processing workflow. In this way, computing the height uncertainty and assuming it as prevalent on others, is can be used to estimate the volume uncertainty.

For an evaluation of the accuracy achieved of the altimetric models produced by means of the photogrammetric process, a comparison was carried out with CloudCompare using the same above-mentioned modalities, between the point clouds produced with the photogrammetric process and the one produced by the laser scanner survey. For this comparison, only some sample areas were selected, characterized by good coverage and geometry of the laser-scanner point cloud.

In the comparison between the photogrammetric point cloud of the vertical UAV acquisition (before called A) and that of the laser scanner, 450,000 points were calculated and compared, obtaining an average value of the differences in Z equal to approximately −1.00 cm and a standard deviation of about 1.36 cm.

The comparison carried out between the photogrammetric point cloud relating to the vertical + oblique flight (C) and the laser scanner acquisition, again on about 450,000 points, provided an average value of the differences in Z equal to about 0.3 cm with a standard deviation of about 1.45 cm.

The results demonstrate that the two scenarios are substantially equivalent; however, compared to practically identical standard deviations, the average value of the deviation is much lower for the point cloud obtained with the vertical + oblique acquisition (C), which in some ways is an indication of a minor systematic component (Figure 27).

On the basis of the above considerations, the evaluation of the accuracy in the computation of the volumes was carried out with reference to the computed values for the vertical + oblique UAV acquisition (C), without prejudice to the substantial equivalence of the results obtained with only the vertical shoot.

Considering that the standard deviation at height ${\sigma }_h$ can be computed by the standard deviations of the height precision of ground points and surface points, for the value to be assigned, it has been considered that both the measurement of the base areas and the stockpiles surfaces were affected by the same biases. The previously determined deviation value was therefore considered for both the base and the surface and the two standard deviations were statistically summed by means of the law of error propagation, thus determining a value of ${\sigma }_h$ equal to ±2.05 cm.

Trials and studies carried out show that the link between the standard deviation of the computed volume and that of the surface model used can be represented by the following equation [39]:

$$\begin{gathered} {\sigma }_V=S · {\sigma }_h \\ \mathrm{S}=\mathrm{area} \end{gathered}$$

(4)

In view of the simplification of considering the calculated standard deviation value for the whole area and applying Equation (4) to the volume values calculated with ArcGIS for the vertical + oblique acquisition scenario (C, see Table 6), the estimated error values have been calculated as indicated in

Table 8. Computed volumes with estimates of the errors.

	V (m3)	S (m2)	${\mathit{\sigma }}_{\mathit{V}}(\pm {\mathit{m}}^3)$	$\frac{{\mathit{\sigma }}_{\mathit{V}}}{\mathit{V}}(\%)$
0	398.40	197	4.04	1.01
1	465.05	221	4.53	0.97
2	391.09	185	3.79	0.97
3	1256.03	457	9.37	0.75
4	74.75	86	1.76	2.36
5	619.57	267	5.47	0.88
6	5559.63	1401	28.72	0.52
7	2060.83	461	9.45	0.46
8	1535.67	753	15.44	1.01

.

The percentage values are between 0.46% and 2.36%, with the maximum value for the volume “4”, and with the exception of this, we have values all below or very close to 1%; this complies what previously stated in Section 2.2.4 about relative errors that are generally more influential for piles with a high surface/volume ratio (see also Figure 9). It must be also considered that the uncertainty in measurement of the ground area has been neglected, as all computations are referred to the same areas, but this can cause significant errors in piles with a high S/V ratio.

In reference to volume “6”, in Table 8, marked in red, the estimated error value is not considered reliable in virtue of the excessive deviation in the model (as already highlighted in Section 5.1). For this reason, the assumption that the standard deviation calculated in the sample is applicable to the entire model is not deemed valid. The calculated value has been presented in any case for demonstration and completeness purposes only.

The results are quite satisfactory and perfectly in line—if not better—than those found in similar experiences and studies; this leads to confirm the perfect applicability, in terms of accuracy, of UAV photogrammetry for computing the volume of bulk materials. Similar experiences in the methods, but carried out on stockpiles deriving from extraction activities, indicate as 3% the value of admissible accuracy established by different regulations pertaining to mining activities [39].

5.3. Evaluation of the Timing

The timing of the activities described, starting from the phase of the photogrammetric survey up to the final determination of the volumes is summarized in

Table 9. Summary of the work times with the UAV photogrammetric survey method.

Work Phase	Times
UAV photogrammetric survey, including the preparation and control phases	2 h (one shoot) 3 h (two shoots)
Pre-processing, alignment, GCPs insertion and collimation, optimisation	5 h (approx. 3 for the computation)
Dense cloud processing	From 4 to 24 h (depending on the hardware)
Editing of dense cloud, creation of mesh + texture, creation of DEM	1 h
Acquisition and exporting delimitation of volumes	1 h
Ground classification and creation of DEM base	1 h
Photoscan: volume computation	30 min
ArcGIS: data importation and volume computation (builder model)	30 min
TOTAL SURVEY AND COMPUTATION PHASE (with the best hardware configuration)	3 h survey; 5 h interactive processing; 7 h automatic processing. 30 min computation

.

The phases pertaining to the survey and the processing of the laser acquisitions have been excluded, considering only the workflow for a periodic monitoring carried out with the described photogrammetric survey. For the same reason, the positioning and survey phases of the GCP coordinates were not included because, as already stated, in case of periodic monitoring, targets to be used as GCPs would be permanently positioned and measured only during the first system set-up or if they have to be repositioned. On this premise, it is reasonable to assume that the stockpiles volume can be computed within 24–36 h after the survey.

6. Conclusions

With the UAV photogrammetry, the volume computation based on the three-dimensional modelling of the surveyed surfaces reached accuracies perfectly adequate for the practical requirements of the case study. The comparison with surface models generated with a laser-scanner survey, which through the extraction of GCPs also provided the absolute reference system, showed average errors of a few centimetres from the value considered as ground truth. These errors in the surface model translate into average percentage errors of approximately 1% on the computed volumes.

Major errors localized in specific areas are due to a non-uniform distribution of the temporary GCPs, despite being adequate for the purposes of the methodological test. In a real implementation of the tested method, the best solution with reference to the GCPs would seem to use fixed targets uniformly distributed in the area and along the boundary, to be measured by other survey methods.

As far as the configuration of the acquisition is concerned, the advantages of an additional oblique camera acquisition do not appear evident; the slight differences in accuracy, even if in favour of the vertical + oblique acquisition, do not appear statistically significant in the face of higher costs and longer survey times.

An element to which the utmost attention must be paid is that of defining the surface of the reference base obviously mostly covered; the solution tested for its modelling, consisting of a classification of the “ground” from the cloud of photogrammetric points, also appears satisfactory and could be improved if the same method is applied to a survey carried out in periods of minimum material accumulation. However, the best solution seems to be that of performing, always in periods of lower accumulation, a high-precision survey aimed at the acquisition with the greatest possible accuracy of all the “fixed” parts of the area in question; this simultaneously provides both the fixed and stable reference of the base, as well as the definition of all the boundary structures which, due to being partially covered in the study case, could not be completely excluded from the computation.