Automatic quality estimation in blending using a 3D stockpile management model

Automatic quality estimation in blending using a 3D stockpile management model

Advanced Engineering Informatics xxx (2015) xxx–xxx Contents lists available at ScienceDirect Advanced Engineering Informatics journal homepage: www...

7MB Sizes 0 Downloads 32 Views

Advanced Engineering Informatics xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Advanced Engineering Informatics journal homepage: www.elsevier.com/locate/aei

Automatic quality estimation in blending using a 3D stockpile management model Shi Zhao a, Tien-Fu Lu a,⇑, Ben Koch b, Alan Hurdsman b a b

School of Mechanical Engineering, The University of Adelaide, SA 5005, Australia Eka Australia, Level 5, 185 Victoria Square, Adelaide, SA 5000, Australia

a r t i c l e

i n f o

Article history: Received 10 February 2015 Received in revised form 10 June 2015 Accepted 12 July 2015 Available online xxxx Keywords: Blending 3D stockpile voxelization BWR automation Mining planning Quality management

a b s t r a c t Stockpile blending is always considered to be an effective method of controlling the quality and maintaining the grade consistency of delivered bulk material, such as iron ore. However, major challenges remain to predict the quality of a stockpile during blending (stacking and reclaiming) operations, because the chemical composition of the ore body is not always available during blending. Consequently, the performance of current stockpile management systems is relatively poorer than expectations. This paper details an innovative modelling approach to estimate the quality of a stockpile during the blending process. The geometric model created from laser scanning data is capable of recording the dynamic shapes of the stockpile using mathematical equations. Therefore, the quality of the stockpile is calculated with a high degree of accuracy when the chemical analysis is completed. Thus, a quality embedded geometric model is created. Furthermore, this geometric model is associated with the kinematic model of a Bucket Wheel Reclaimer (BWR) to achieve precise auto-control in reclaiming operations. The link between these two models also supports the quality of every single cut accomplished by the BWR, to be calculated in advance. Using the calculation results in conjunction with precise machine control, the output grade can be predicted, planned and controlled. This will optimise the efficiency and effectiveness of stockpile blending. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction A major concern in iron ore exporting is the delivery of a relatively constant quality that is close to the target grade. Mining practices, ore processing, railing, blending and ship loading operations construct a multi-level and multi-parameter control system for such a purpose. Out of these operations, stockpile blending is widely accepted as an effective method to adjust the quality and improve the homogeneity of the final product during the handling process. These two functions are achieved through stacking a stockpile into different geometric shapes layer-by-layer and then reclaiming the stockpile slice-by-slice. Thus, materials in separate layers which have different chemical compositions are mixed and the quality fluctuations are reduced. Stockpile management systems are designed to assist the decision making in blending by researchers from different disciplines. Generally, these models can be classified into two types:

⇑ Corresponding author. Tel.: +61 8 8303 3556; fax: +61 8 8303 4367. E-mail addresses: [email protected] (S. Zhao), [email protected] (T.-F. Lu), [email protected] (B. Koch), [email protected] (A. Hurdsman).

mathematical models [1–3] and geometric models [4–6]. A more detailed review of these models can be found in the authors’ previous study [7]. A common idea in these modelling approaches is to stack newly arrived ore selectively, each according to its quality composition, so as to optimise the homogeneity in the quality of that stockpile. However, an intrinsic challenge in such ‘selective stacking’ approaches is that the chemical properties of the ore may be not available beforehand due to the limitations of current sampling techniques. There is always a delay of several hours in obtaining the most accurate analytical results according to current sampling standards. As a result, either received iron ore is held until the quality data are available, which is less likely to happen because it will lower the system efficiency and occupy the storage space, or it is stacked onto a stockpile using established quality information. However, such information may be inaccurate and out-of-date. For instance, when the train delivery of iron ore is crushed and screened into lump and fines products, the separation will cause systematic differences in the chemical compositions of the two products. Both pieces of quality information are needed to direct these two products to specific stockpiles. Conversely, there is only one set of data which is obtained before crushing and screening. The use of such inaccurate quality data may cause

http://dx.doi.org/10.1016/j.aei.2015.07.002 1474-0346/Ó 2015 Elsevier Ltd. All rights reserved.

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

2

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

the optimization algorithm to fail and the quality of a stockpile not to be known exactly after stacking. Another chance to adjust the quality of delivered iron ore occurs at the reclaiming phase. Reclaiming operations are conducted at different sections of a stockpile or across stockpiles in order to adjust recovered material according to the final objective. However, due to the lack of accurate quality knowledge of the stockpile(s), the operation has to be suspended frequently while waiting for the chemical analysis results. If any changes in quality are needed, the BWR(s) has to move to another position to recover different grade material. Again, the decision is made through imperfect quality information of the stockpile. Furthermore, the reclaiming has to be activated manually by an experienced operator to avoid any collision between the stockpile and the wheel of the BWR. Thus, such semi-automatic control inevitably introduces more uncertainties into the quality adjustment process. Consequently, the adjustment may not balance the bias in quality as expected. Although BWRs are operated semi-automatically and the reclaiming efficiency is relatively poorer than expectations, it is surprising that only a few studies have reviewed this area. To land the bucket wheel (BW) onto a stockpile automatically, Choi et al. extracted a 2D contour map from 3D laser scanning data and evaluated the points on the 2D contour lines over the slewing trajectory of the BW [8,9]. The validated landing points will not cause a collision between the wheel body and the contour line of the stockpile throughout the reclaiming process. The joint angles of the BWR are calculated from the landing point using inverse kinematic equations. Such idea was adopted by Lee et al. in 2006. They built a contour map from scanning data to optimise the trajectory of a BWR with regard to the safety, energy consumption and transfer time [10]. Clearly, these approaches are not focused on optimising the quality control of the blending operation. Additionally, the boom and the bucket wheel are assumed to be on the same slewing plane in these papers. Thus, the 3D collision detection problem is simplified as a 2D case. Conversely, in practice, to reduce the cutting resistance, the bucket wheel is tilted at a small angle. Therefore, such 2D approximation may be not able to identify all potential collision points in 3D space. To improve the efficiency of the quality control in reclaiming, Lu and Myo described an isosceles triangular prism stockpile model for BWR trajectory optimization. They partitioned the model into a number of sections (triangular/trapezoidal prisms) which are called voxels in their paper and assumed the quality variations in these voxels satisfied a normal distribution [11,12]. Together with other operational constrains, they proposed a trajectory optimization algorithm that achieved the customer-required quality with minimum energy consumption. However, the success of their approach is based on the assumption that the chemical composition in each voxel space is known exactly. This assumption is not quite appropriate regarding the issue existing in the previously discussed stacking optimisation models. Meanwhile, the triangular prism model is far from a good representation of a real stockpile and the shapes of these voxels are quite different from the crescent-shaped cuts generated by a BWR. There is definitely a lack of a novel stockpile management model that is able to link the stockpile and BWR models together to calculate the quality during the blending process in real-time and guide the BWR operations. To summarise, current stockpile models are insufficient for accurate quality calculations and efficient quality control. The main reasons for such a situation include: the geometric shapes of a stockpile are not acquirable in real-time. The stockpile models are isolated from the BWR model in terms of quality control. Therefore, this paper presents a novel geometric model that is embedded with quality information. Through associating this integrated model with a BWR model, it is possible to control a BWR

automatically and estimate the quality of the reclaimed material for each cut. As a result, it is possible to predict the tonnage and quality grade with a great degree of accuracy and continually adjust the reclaiming according to the target grade instead of using current selective stacking approaches. In other words, it is possible to produce the desired product quality through controlling the BWR movement proactively during reclaiming operations. The use of this model will change current stockpile management fundamentally and improve bulk material handling efficiency significantly. The rest of this paper is organised as follows: Section 2 describes the proposed BWR and stockpile models respectively. Section 3 details the algorithm that integrates the quality information into the geometric model and how to achieve automatic operation of a BWR using the 3D stockpile management model. Section 4 presents two case studies regarding two groups of datasets. Section 5 concludes the paper, together with a discussion of future work. 2. 3D stockpile management system modelling The proposed 3D stockpile management model contains two components: a BWR model and a 3D stockpile model. This section details the modelling algorithms. 2.1. BWR modelling A boom-type Bucket Wheel Reclaimer (BWR), which combines the features of a stacker and a reclaimer, has been widely used for stacking/reclaiming in iron ore handling. In stacking mode, it transfers the material through the conveyor belt along the boom and allows the material to fall at the end of the boom. In reclaiming mode, the boom slews across the body of a stockpile and the buckets on the rotating wheel scoop materials from the pile simultaneously. A kinematic model deduced by Lu [13] was used to represent the movement of a BWR in 3D space. This model has been proven to be accurate and efficient for both initial and final position calculations [14]. Using this model, the position of the end-effector (the centre of the bucket wheel) can be calculated through the joints’ angles, and vice versa. Assuming that the BW rotates at a constant speed and the material is being scooped evenly by the buckets during the cutting phase, it is possible to simplify the BW as a circular ring. Meanwhile, if the vibration of the boom during the slewing motion is ignored, the BW can be considered to rotate perfectly according to its lateral axis. Under these assumptions, if the rotation axis of the BW is perpendicular to the boom surface, the slewing motion of the boom will result in a non-degenerate torus surface out of the circular ring. The radius of the tube is equal to the radius of the BW and the distance to the axis of revolution can be calculated using the kinematic equations. If the BW is titled at an angle of b, the projection of the circular ring to its rotating plane results in an ellipse and the revolving of the BW will generate an elliptic torus. The horizontal and vertical semi-axis of the elliptical cross-section is determined by the tile angle and the radius of the BW. Fig. 1 plots a section of two such types of torus, separately. 2.2. Stockpile modelling To build an accurate and up-to-date 3D model, a mobile device is essential to allow the sensor to be mounted on and measure the stockpile’s surface as it is stacked/reclaimed. According to the literature survey conducted by the authors in the previous paper [15], a scanning system developed by a company named QMASTOR is built through mounting a 3D Lidar at the back of a vehicle [16].

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

3

Fig. 1. Ideal surface generated through superimposing the rotation motion of the BW and the slewing motion of the boom. The magenta line with green squares indicates the boom of the BWR. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

To measure the stockpile, the vehicle has to be driven around the stockyard. Additionally, a company called Indurad employs a radar sensor to measure the profile of the stockpile. An extra supporting frame is constructed on top of the stockpile to drive the radar [17]. Due to the protection of intellectual properties, the authors have no access to the detailed design of either system. The above knowledge is acquired from their websites, images, videos and contact emails. However, it is apparent that these two systems require extra hardware and manpower for stockpile scanning missions. The solution proposed by the authors is to upgrade a BWR into a mobile scanning device. When a 2D laser scanner is mounted at the end of the boom of a BWR and faced toward a stockpile, the motion of the boom is exploited to move the scanning plane and produce a de facto 3D scanner that measures the profiles of a stockpile simultaneously during stacking and reclaiming operations. This is a cost-effective solution because a 2D laser scanner is much cheaper than a 3D Lidar. Also, such a scanning mechanism is more flexible and efficient than current commercial solutions that have been identified in the literature. Currently, most BWRs operated in Australia have on-board encoder systems to measure the position of the buck wheel for the purpose of automation and collision avoidance. However, encoders’ measurements suffer from quantization errors and degrade without limitation over time. Errors by travelling encoder(s) are eliminated by placing many calibration points along the BWR’s moving track. However, no calibration can be applied to the slewing and luffing angular encoders, due to the technical difficulties in mounting extra calibration devices to such large-scale machinery. As a consequence, the estimated position errors of the end-effector are more than 30 cm, resulting in a significant corruption of the Lidar measurements. Furthermore, this error range does not contain the sources caused by the imperfect mechanical linkage and the vibration of the boom during machine movements. Such positioning accuracy is not favourable for stockpile modelling. To map the environment using a mobile robot, the mobile robot itself needs to be localised in advance. Localization denotes the capability of a robot to elaborate the sensor measures and find out its position with respect to a coordinate system. Localization is a fundamental problem in mobile robots. Generally, Kalman Filter (KF) based sensor data fusion techniques are widely used for outdoor mobile robot localization. For this application, the authors proposed an Unscented Kalman Filter (UKF) based algorithm to fuse the data provided by on-board encoders and an extra GPS mounting, together with the 2D laser scanner (see Fig. 2). In contrast with encoders, the accuracy of the GPS neither degrades over time nor distance. Due to their complementary characters, the integration of GPS and encoders through Kalman filtering techniques will be a good match to improve the position estimations. A UKF can be considered as a KF based on an unscented transformation (UT). The UT generates a set of weighted sigma points,

Fig. 2. Mobile laser scanning device concept design. The encoders installed to measure the travelling, slewing and luffing motions are highlighted in green. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

which are located at the mean and symmetrically along the main axes of the covariance, to represent the Gaussian distribution and linearize the nonlinear system function. The use of these weighted sigma points in the UT propagates more accurate mean and covariance results than the use of the first-order Taylor expansion in linearization of an Extended Kalman Filter (EKF). Readers who are interested in the implementation of UKF based sensor data fusions for BWR localization can find a detailed description in the author’s previous paper [18]. According to the simulation results, the overall positioning accuracy is better than 20 cm when the localization algorithm is applied to a BWR with 10 m in height and 50 m in boom length. After consulting with the industry partner, MatrixGroup, the authors were informed that it is currently not possible to mount a laser scanner on a BWR and scan a real stockpile because it is hard to find the right experimental venue and a spare BWR for the research. Therefore, a 3DOF laser scanning system was built to scan a scaled-down stockpile model in the laboratory environment [19]. The LMS200 laser range finder manufactured by Sick was mounted on a pan and title platform (QPT-50). The platform was fixed on a roller which was driven by a motor through the roller chain linearly. Point cloud data acquired from the scanning system are organised into a three-column matrix. Each column represents an axis in the Cartesian coordinate system (PLaser (pLaser_x, pLaser_y, pLaser_z)). The letter after the underscore in the subscript means the axis of abscissas. Fig. 3 illustrates the CAD drawing of the 3DOF laser scanner used in the laboratory environment. The scanning mechanism of this indoor scanner is similar to the proposed mobile laser scanner.

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

4

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

knot vector U = {u0, . . . , um}. Each ui is referred to as a knot and ui 6 ui+1, i = 0, . . . , m  1. The recursive B-Spline function is called Cox-De Boor relations [21]. A spline in the vector u is:

 Ni;0 ðuÞ ¼

Ni;di ðuÞ ¼

Fig. 3. A 3DOF laser scanning system designed for laboratory use.

The surface modelling aims to reconstruct the original shape of the stockpile using mathematical equations. Pre-processing is applied to the raw data to segment the stockpile region from the background and filter out the measurement noise before modelling. Image processing techniques are used to extract the stockpile boundaries from the point cloud [7]. Once the boundaries are detected, points outside the stockpile region are eliminated. Meanwhile, based on the boundary detection results, the point cloud data in the scanning reference frame can be converted into its local coordinates, which represent real elevations of the stockpile to the ground surface. A wireframe model presented in the previous work [20] is used to smooth and store the measurement data. In this model, a profile of the stockpile along the scanning direction is approximated by a 4th order Fourier series function. The main reason of choosing Fourier descriptors is that they smooth the data without obvious shrinkage effects. Fig. 4 plots the point cloud data before and after the pre-processing. A 3  3 Gaussian filter can be applied to smooth the data further in the z-axis before creating a surface model. However, this filter is optional. The reason will be discussed in the results section. A B-Spline tensor product surface S(u, v), defined by a 2D control net Qi,j with 0 6 i 6 m and 0 6 j 6 n, is a vector-valued function of two parameters, u and v (in this case, they represent the x and y axis respectively): m X n X Sðu; v Þ ¼ Ni;di ðuÞNj;dj ðv ÞQ i;j

ð1Þ

i¼0 j¼0

where the di and dj are the degrees for the surface function and 1 6 di 6 m, 1 6 dj 6 n. The function N i;di ðuÞ and N j;dj ðv Þ are the dith and djth-degree B-Spline basis functions. They are defined reclusively through a sequence of non-decreasing real numbers called

1; if ui 6 u < uiþ1 0;

ð2Þ

otherwise

uiþdi þ1  u u  ui Ni;di 1 ðuÞ þ Niþ1;di 1 ðuÞ uiþdi  ui uiþdi þ1  uiþ1

ð3Þ

N i;0 ðuÞ is a step function, equal to zero everywhere except on the half open interval u 2 ½ui ; uiþ1 Þ, and for a di (di > 0), Ni;d ðuÞ is a linear combination of two (di  1)-degree basis functions [22]. According to such a definition, the knot vector U determines the functions Ni;di ðuÞ completely once the degree of the function is fixed. Similar definitions can be also extended to the vector v. The knot vector can be either periodic or non-periodic (open). In terms of the geometric shape of a real stockpile, the knot vector used for the surface model is non-periodic and is determined by the centripetal method [23]. The idea of the centripetal method is that the normal acceleration (i.e., centripetal force) should not be too large and should be proportional to the change in angle when a car is driven through a sharp corner. Supposing there are l measurements (P0, P1, . . . , Pl) along the x-axis of the stockpile region [a, b] that obtained from previous boundary detection results and a = min(Pl_x), b = max(Pl_x), the distance between two adjacent data points is measured by |Dk  Dk1|a rather than the true distance |Dk  Dk1| (usually, a = 1/2 for the square root). The total length of the data polygon L is:



l X jDi  Di1 j1=2

ð4Þ

i¼1

The ratio of the distance from the data point D0 to the point Dk, denoted as Lk (k = 1, . . . , m + d  1) over the total length, is:

Pk

i¼1 jDi

Lk ¼

 Di1 j1=2 L

ð5Þ

Therefore, the knots for U are:

u0 ¼ a;

uk ¼ a þ Lk ða  bÞ;

um ¼ b

Both the coefficients of the spline surface (the control point Qi,j) and the point cloud data that represented the stockpile region (P) can be arranged in a matrix format, respectively:

2

Q 00 6 . Q ¼6 4 .. Q m0

3 Q 0n .. 7 7 . 5;    Q mn  .. .

2

P00 6 . P¼6 4 .. Q p0

3 P 0q .. 7 7 . 5    Q pq  .. .

The Q is an (m + 1)  (n + 1) matrix, while P is a (p + 1)  (q + 1) matrix. The least square error function that measures the square distance between the B-spline surface and laser scanning points is:

Fig. 4. A point cloud of a stockpile after the pre-processing. The red line shows the boundary of the stockpile. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

5

Fig. 5. A surface patch rendered from the point cloud using a tensor product spline. Profiles along the longitudinal direction, generated from the Fourier model, are plotted in grey dots.

Fig. 6. Landing the bucket wheel on the stockpile. (a) The maximum cutting depth (H) of the BW in regard to the bank surface. u is the angle at which material starts to fall from the bucket and the b is the titled angle of the bucket wheel. The shadowed area indicates the material that is recovered by the buckets. (b) An ideal case when determining the BW centre from the surface normal. The normal vector Ns of the least square surface at the point pq is tangent vector U. U is in the opposite direction of the slewing direction of the bucket wheel at pq. Vector V is directed towards the slewing axis of the BWR. The angle between the Ns and U is a and a = 90°.

2  p q  m X n 1 X X X  ^ EðQ Þ ¼ N i;di ðuk0 ÞNj;dj ðv k1 ÞQ ij  Pk0 k1    2 k ¼0k ¼0 i¼0 j¼0 0

ð6Þ

1

The solution of (6) is the control points (coefficients) that minimize b Þ. The global minimum occurs when all its first order partial Eð Q derivatives are zero: ! p X q m X n X X @E ¼ Ni;di ðuk0 ÞNj;dj ðv k1 ÞQ ij  P k0 k1 Nk0 ;di ðuk0 ÞNk1 ;dj ðv k1 Þ @Q k0;k1 k ¼0k ¼0 i¼0 j¼0 0

¼

1

p X q X m X n X Ni;di ðuk0 ÞNj;dj ðv k1 ÞNk0 ;di ðuk0 ÞNk1 ;dj ðv k1 ÞQ ij k0 ¼0k1 ¼0 i¼0 j¼0



p X q X

Nk0 ;di ðuk0 ÞNk1 ;dj ðv k1 Þ  P k0 k1

k0 ¼0k1 ¼0

where ak0 ;i ¼ N k0 ;di ðuk0 Þ and ak1 ;j ¼ N k1 ;dj ðuk1 Þ, a matrix format of this equation is:

@E ¼ AT AQBT B  AT PB ¼ 0 @Q

ð7Þ

where A is a (p + 1)  (m + 1) and B is a (q + 1)  (n + 1) matrix, respectively.The coefficients matrix Q is: 1

Q ¼ ½ðAT AÞ AT P½BðBT BÞBT 

T

ð8Þ

Two sets of basic functions in the u and v directions are created after fitting. When needed, a desired point data that represents the stockpile surface can be generated through the B-spline surface function (see Fig. 5). 3. 3D stockpile management model Because both the stockpile and BWR models have their own reference frames, a new global coordinate system (right-handed) is defined to integrate the two. The centre of the new coordinate is located at the centre of the BWR kinematic model. The positive z-axis, normal to the crawler level, forms the slewing axis. The positive y-axis is in the direction of the BWR’s advancement. The positive x-axis is normal to the y–z plane. Namely, once the distance between the stockpile boundary and the BWR centre is known, a 3D stockpile management model is generated after the coordinate transformation. The following section details the two most important applications for this model: machine automation and quality estimation in reclaiming. 3.1. BWR automation When a stockpile is recovered using the bench reclaiming (or terrace cutting) method, an operator selects a landing point on

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

6

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

Fig. 7. Evaluating the landing point set in the stockpile management model. The stockpile is partitioned into a 19  10 grid. Each cell in the square grid is 80 mm. The query points are highlighted in green. The effective cutting region, determined by the maximum cutting height, is enclosed by the magenta curve. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8. Updating the geometric shape of the stockpile model. The intersecting curve between the torus and the stockpile highlights an enclosed region. In this region, points on the stockpile’s surface will be replaced by the points on the torus.

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

7

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

the stockpile surface manually and approaches the BW to that point with great care. When the contact is about to happen, the operator activates the circular motion of the buckets to make sure the buckets can pick up materials from the stockpile simultaneously with the landing. Normally, the rest of the operations are controlled automatically until the height of the BW is changed. The slewing range is determined by the feedback from the loading sensor. The lateral motion begins to stop until the resistance force is no longer detected by the sensor. After a single cut, the BWR is advanced with a fixed cutting depth (i.e. 0.9  bucket width) to conduct a new cut. The BW is lowered to start a new reclaim when all the materials at that terrace are recovered. Again, the operator needs to choose the correct landing point for the new reclaim. Therefore, BWR automation contains two tasks: landing the BW onto the stockpile and predicating the slewing range for each cut. More importantly, landing is crucial for future quality control because the cutting depth of the BW is fixed in most cases. The problem of BWR automation in this paper is stated as: given a set of candidate points Pc on the stockpile surface, validate a subset Pv that must not cause collision between the wheel body and the pile of all the cuts in one reclaim. Then, calculate the initial joint angles of the BWR (d1, h2, h3) from the subset Pv for landing the BW automatically onto the stockpile. Additionally, calculate the slewing range of each cut from Pv required for the BWR operation control. A cut means a revolving motion of the boom and a reclaim is a collection of cuts with an equal cutting height.

3.1.1. Landing point estimation The stockpile model contains an infinite number of points and every point could be a candidate point. It is neither feasible nor necessary to assess all these points. In this paper, the candidate point set Pc is determined by three thresholds: the effective cutting depth, the facing direction and the square grid spacing. The effective cutting depth (De) is defined based on the maximum cutting height of a BWR, as shown in Fig. 6a. Assuming that material inside the bucket starts to fall due to the gravity at a certain angle u, the maximum cutting height H approximately equals (Rv + Rksin u)cos b and De = max(PL_z)  H. Additionally, because the BWR is normally mounted on rails or tracks, only one side of the stockpile is available for the landing mission. The first two thresholds can be summarized as: only those points within the effective cutting region and facing the BWR mounting direction will be considered for landing. However, there are still infinite points after thresholding. To obtain a finite set, the stockpile region on the xy plane is partitioned into a group of contiguous square cells and the centres of these grids are superimposed on the stockpile surface. Thus, a finite set Pc is created after such a grid reference method. Obviously, a small grid spacing (high resolution) will generate more points in the set. To calculate the joint positions of the BWR for a given query point pq in Pc using inverse kinematic equations, the position of the end-effector needs to be identified. Due to the rotation of the BW, the exact contact point between the BW and the stockpile surface cannot be solved using the kinematics of the BWR directly. However, because the purpose of the auto-landing is not to locate the exact point on the BW that is tangent to the stockpile surface at pq, instead of adding more constraint functions into the model, a simple way to approximate the position of the BW centre Cbw is by assuming that the vector between two points (Cbw and pq) is in accordance with the surface normal at the point pq (Fig. 6b). Therefore, Cbw = pq + Rvns, while the normal vector ns at pq is calculated using the stockpile model. Once the Cbw is known, the d1, h2, h3 can be solved using inverse kinematic equations. The definition for d1, h2, h3 and L2  L4 can be found in Lu’s paper [14] and are shown in Fig. 7.

To detect the collision between the wheel and the stockpile in one cut is to detect whether the points on the stockpile surface are also inside the torus generated by the revolving motion of the BW. If the cross-section of the toroidal surface is a circle, a point ps (ps_x, ps_y, ps_z) on the stockpile surface is inside the toroidal surface if:

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 Rm  ðps x  x0 Þ2 þ ðps y  y0 Þ2 þ ðps z  z0 Þ2 < R2k

ð9Þ

where Rm is the major radius of the torus and Rm = L5sin h3 + L4 and Rk is the radius of the wheel which is equal to the minor radius of the torus. x0, y0 and z0 is the centre of the torus and is calculated from the stockpile management model. When the cross-section is an ellipse, a point is inside the elliptic torus if:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðps x  x0 Þ2 þ ðps y  y0 Þ2  Rm Rk cos b

þ

ðps z  z0 Þ2 <1 Rk

ð10Þ

where b is the inclined angle of the BW. 3.1.2. Slewing range estimation and model updating If no collision is detected in a single cut, the slewing angle of the BWR will be calculated and the recovered region needs to be replaced by the helical-shaped cutting surface (a part of the torus) for the next assessment. For these two purposes, the stockpile is presented as a set of points in 3D space calculated from its spline model. The elimination of recovered points from stockpile is quite easy because Eqs. (9) and (10) are still validated. The only change is to replace Rk by Rv. The cutting surface is the region that is enclosed by the intersecting curve between the torus (T(s, t)) and the stockpile (S(u, v)), supposing that the torus is defined topologically as the product of two circles s and t. The intersecting curve is an algebraic equation T(s, t)  S(u, v) = 0. Once the equation is solved, the cutting surface can be obtained to update the stockpile model. However, the relatively high degree of the algebraic representation may lead to an algebraic curve of even higher degree. Solving such a high-degree algebraic curve is a complex problem. For this reason, a numerical approach is proposed to determine the intersecting curve and the slewing range. This numerical approach separates the continuous slewing motion into a group of discrete steps. At each step, the centre of the BW is calculated using kinematic equations and represented as a set of points Pbw. These points define a circle (b = 0) or an ellipse (b – 0) in 3D space. Theoretically, the intersecting point(s) q (q 2 Pbw ) between the circle/ellipse and stockpile surface is S(q_x, q_y)  q_z = 0. In other words, the Euclidean distance D between point(s) q and the stockpile surface S(u, v) is zero. Again, as the algebraic degree and genus of this equation is too high to be solved mathematically, a numerical approximation would be finding out the points from Pbw that have a minimum Euclidean distance to the stockpile surface and the minimum distance min(D)  0. If min(D) is much larger than zero, the circle/ellipse is believed not to be intersected with the stockpile. Because the BW will only intersect with the stockpile at two points at most, this numerical approach is sufficient to find out the right solutions. Once the intersecting points are detected, Pbw can be separated into two subsets: inside and outside the stockpile. Points inside the stockpile and the intersecting points will be inserted into the stockpile model. Points outside the stockpile are eliminated. After searching all the slewing steps, the actual angles that the BW cuts into and moves out of the stockpile can be obtained. This slewing range will be used for the BWR’s auto-control. Fig. 8 plots the shape of the stockpile after the first two cuts. The elimination

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

8

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

Fig. 9. Calculating the volume of a voxel (cut) that contains one QVO, using the QMC method.

Fig. 10. Partitioning a voxel into octants to simulate the quality variation.

and model updating procedure proceeds iteratively until a reclaim is completed or a collision has been detected.

3.2. Quality estimation in blending The quality of a particular element acquired after the chemical analysis is the percentage by weight with the mean and standard deviation of the entire compound. Given the truth that the sizes of ore particles of the same product are generally in the same range, the density of the iron ore can be considered as a constant. Therefore, the quality of a certain chemical component in a 3D space is determined by the volume of that 3D space. As regards a stockpile that contains multiple layers, the quality of the component is linked with the volumes of these layers because the chemical compositions of these layers are different. Since each layer is scanned by the laser and represented by a spline surface function, the volume under the surface can be calculated using a double or triple integral. Given that the calculation procedure is available in any calculus book, the authors will not give any further descriptions of the process in this paper. Similarly, the quality of the reclaimed material is determined by both the volume of the cut and the chemical composition of the stockpile layers inside the cut. The calculation algorithm proposed here aims to link these two correlated factors using geometric information. In this paper, a sickle-shaped cut caused by the BWR is considered as a quality unit, which is called a voxel. A voxel may contain materials from a number of layers and the chemical composition may vary from layer to layer. Consequently, the volume of each portion of the layer needs to be calculated separately. If the intersection between a voxel and a layer is called a quality volume object (QVO), a voxel intersecting with multiple layers will have multiple QVOs. Since each QVO belongs to a layer, it shares

the same physical properties (i.e. chemical composition and density) with that layer. Due to the variability in shapes of these objects, it is difficult to program a universal algorithm to determine the integral region and perform the calculation automatically. Additionally, even if such a program were developed, the complexity of the code would also prolong the calculation time. For these reasons, a Monte Carlo based integration method is selected for this study. The Monte Carlo integration uses randomly generated points to approximate volume integrals of a 3D object. These points are randomly distributed inside a box enclosing an object of interest. The volume of the object is estimated as:

V obj ¼ V box ðNin =NÞ

ð11Þ

where Vbox is the volume of the 3D bounding box of the object, Nin is the number of points that are contained within the object and N is the number of points generated in a random process. A major source of error for such a method is the random nature of the sampling because the pseudo-random points generated from the computer may not uniformly distribute inside the bounding box. Quasi-Monte Carlo (QMC) techniques replace these pseudo-random numbers by deterministic sequences, which are constructed explicitly to minimize clumping, through reducing the discrepancy of the point set. The use of such low-discrepancy numbers rather than pseudo-random numbers in integral calculations increases the convergence rate and reduces the expected relative error. As pointed out by Davies and Martin [24], the error of the approximation by QMC methods is O(N1log3 N). The implementation of the QMC method is quite simple. After generating a low-discrepancy set using Sobol’s method [25], these points are shifted into the bounding box and then Nin is determined through point-membership classification based on the stockpile’s surface and torus functions (1), (9) or (10).

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

9

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

Fig. 11. Quality calculation using QMC method for the 2nd cut. The 1st object is in dark red and the 2nd object is in light cyan. The octant intersecting with two stockpile layers is highlighted in yellow. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The implementation of the QMC method is quite straightforward, as shown in Fig. 9. Given a voxel whose upper surface is a part of the stockpile and lower surface is a part of the torus, a bounding box is created first, based on its geometric shape. Then, a low-discrepancy set with N points, generated by Sobol’s method, is shifted into the bounding box using geometric transformations. Subsequently, based on algebraic functions of the stockpile and torus, the points are separated and the Nin is obtained. Finally, the volume of the voxel Vvoxel is calculated using (11). If a voxel only contains one QVO, and both the density of the QVO (q) and the weight, expressed as a percentage of a certain element (w), are known, the quality of such an element in the voxel Qvoxel is equal to qVvoxelw. Because the quality of iron ore is evaluated by a number of elements, W is constructed as a row vector with we elements. For example, in this paper, the quality of iron ore is evaluated by three elements: Fe, Si and Al. W is a 1  3 vector [wFe, wSi, wAl]. At this stage, it is assumed that the QVO (or the layer) is homogenous in quality (nonhomogeneous cases will be discussed later). If a voxel contains vq QVOs, such calculations will be performed vq times to obtain the quality of the voxel. In turn, a vq  we weight matrix is constructed and the quality of such a voxel is:

Qc ¼

M X

qi V box

i¼1

  NinðiÞ W ði;j¼1:3Þ N

ð12Þ

where qj is the density of the jth QVO and W(j,:) means the jth row of the weight matrix. It is necessary to reflect on the quality fluctuations in the chemical composition during the calculation because a homogeneous layer does not exist in reality. Due to trade secrets, the authors cannot access the real sampling or analytical data. However, with a common understanding that the previous handling processes will improve the homogeneity of the stacked material and the variation in grade will fluctuate slowly compared with the length of the stockpile, the authors assume fluctuations in the quality of a stockpile are normally distributed. Additionally, such normal distribution is a discrete and one-dimensional Gaussian process along

Table 2 BWR parameters for the case study.

Case 1 Case 2

L2

L3

L4

L5

Rk

Rv

u (°)

b

6 cm 6m

5 cm 5m

5 cm 5m

50 cm 50 m

4 cm 7m

5 cm 10 m

20 15

0 5°

the stacking direction. In other words, a stockpile can be divided into a number of segments along the longitude locations (y-axis). Thus, it is possible to assume that the stockpile is partitioned by a sequence of adjacent cuboids and each pair of adjacent faces is overlapped. The quality of such a cuboid is homogeneous and quality variations of these cuboids conform to a normal distribution. To simulate the quality of a stockpile, a subroutine takes the mean and standard deviation of every single element (i.e. Fe, Al or Si) as an input and generates a large amount of random numbers that follow the specified Gaussian distribution. These numbers construct a sampling database and the program will pick up values from the database randomly to represent chemical variations during the computation. To link the geometric information with the chemical variation when calculating the quality of a voxel, the bounding box of the voxel is partitioned into os octants using Octree decomposition [26]. Octants which have the same y-coordinate, are grouped into a segment and the quality composition of a segment is assumed to be homogeneous. In other words, a layer inside a segment should have same quality composition. Meanwhile, the decomposition also subdivides a voxel into os sections and each section is encompassed by an octant (see Fig. 10b). Using a similar grouping strategy as with the octants, sections in the voxel are grouped into og cuboids (segments). Because a voxel is further partitioned, another dimension is added to the weight matrix W. A 3D matrix W [os  vq  we] means the weight matrix has os sections, vq QVOs and we chemical components. The quality of the ith section that contains vq QVOs, calculated using the QMC method, is:

Qj ¼

M X

qi V oct

i¼1

Table 1 Quality data used for the case study. Layer

Density

Fe (%)

Si (%)

Al (%)

Mean

Var.

Mean

Var.

Mean

Var.

Case 1 1 2 3

2.61 g/cm3 2.82 g/cm3 2.75 g/cm3

58.68 57.77 58.19

7.06 6.82 7.13

4.38 3.96 4.54

1.45 1.26 1.32

1.54 2.31 1.66

0.84 0.96 0.75

Case 2 1 2 3 4

5.35 kg/m3 5.35 kg/m3 5.35 kg/m3 5.35 kg/m3

56.23 60.77 51.43 62.05

6.56 4.32 5.18 7.46

4.98 5.33 7.66 5.08

1.53 1.23 1.08 1.34

0.24 0.35 0.41 0.21

0.081 0.063 0.023 0.015

  Ninði;jÞ W i;j;k¼1:3 N

ð13Þ

where Voct is the volume of the octant obtained after octree decomposition. Nin(i,j) means the number of the points within the (ith, jth) QVO and k is the index of the chemical element in the weight matrix. k = 1:3 means the calculation includes all the chemical components. The quality of a voxel that contains os sections is:

Qc ¼

" S M X X j¼1

#   Ninði;jÞ W ði;j;k¼1:3Þ qi V oct N i¼1

ð14Þ

According to previous assumptions in the stockpile quality distribution, octants that have the same y-coordinate form a segment and a layer inside this segment should have the same quality composition. Thus, when the program chooses the quality value from

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

10

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

Fig. 12. A stockpile model created from the point cloud with a 0.25° scanning resolution. (a) Raw data. (b) Range data after the pre-processing. (c) Stockpile surface model. Layers are indexed based on their stacking sequence.

Table 3 Modelling error and execution time using 4th order spline. Layer

Number of points

1 2 3

9704 14,034 17,543

Without filtering

Gaussian smooth filtering

Maximum error (mm)

Average error (mm)

Execution time (s)

Maximum error (mm)

Average error (mm)

Execution time (s)

20.10 12.33 11.26

1.80 0.92 0.70

10.25 12.42 14.24

10.02 6.10 9.47

0.31 0.27 0.30

11.32 13.58 14.86

the sampling database, it is only based on the number of the segments, not the number of octants. Two examples are given in the following section to explain the calculation procedure. Example 1. Estimating the quality of the first cut as shown in Fig. 8. This voxel only contains one QVO and the density of the material q is assumed to be a constant. After octree decomposition, 18 octants are obtained and these octants are classified into 5

segments along the y-coordinate (see Fig. 10). Thus, the voxel is partitioned into 18 sections and the weight matrix W is a (3  18  3) matrix. Taking the first two octants at the bottom left, whose centre is (328, 586, 168) and (348, 586, 168) for example, the quality for each section inside the octant is Q1 = qVoct(N(1,1)/N)W1,1,1:3 and Q2 = qVoct(N(1,2)/N)W1,2,1:3, respectively. Because the 1st and 2nd sections belong to the same segment, W1,1,1:3 is equal to W1,2,1:3. With the 3rd section, the

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

11

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

authors consider this stockpile contains four layers and each layer associates with one point cloud. The quality data used for the case studies are also provided by the MatrixGroup and are listed in Table 1. The BWR models simulated for reclaiming operations are listed in Table 2. All these parameters are, of course, subject to change to reflect the real operational environment. 4.1. Laboratory model

Fig. 13. Residual of the 3rd layer without applying the smooth filter. For display purposes, the authors only highlight the boundaries of the conical ends.

quality Q3 = qVoct(N(2,1)/N)W2,1,1:3. Because there is only one QVO, j always equals 1. The quality of the voxel is the sum of all sections P18 i¼1 Q j .

Example 2. Computing the quality of the second cut, as shown in Fig. 8. The shape of the voxel is shown in Fig. 11 and the boundaries of the cut are highlighted in magenta. This voxel contains materials from two layers. Thus, it has two QVOs (see Fig. 11). The voxel is divided into 16 sections by octants. Regarding its geometric position, a section may have one or two QVOs. If the section only has one QVO object, similar calculation procedures can be used. For example, the quality of the bottom left section is Q1 = q1Voct(N(1,1)/N)W1,1,1:3. If a section contains two QVOs (the one with an octant highlighted in yellow, as shown in Fig. 11), the point-membership classification needs to be executed twice. The quality of the section is: Q1 = q1Voct(N(6,1)/N) W6,1,1:3 + q2Voct(N(6,2)/N)W6,2,1:3. 4. Case study In this section, two case studies relying on two data sources are presented to demonstrate the proposed modelling and quality calculation algorithm. In the first case, the point cloud data are collected from a laboratory scale chevron stockpile which contains three layers scanned by the 3DOF laser scanning system. The scanning range of the laser is 100° and the angular resolution is 0.25°. More detailed experimental settings and procedures are available from the authors’ previous work [20]. Data used for the second case were provided by MatrixGroup, and represent a collection of measurements of a full-scale stockpile. A 70 m  155 m stockyard is sampled by surveyors at a 1 m  1 m quadrat. The shapes of the stockpile after four stacking operations are recorded and the dimensions of the stockpile in the last point cloud are 50 m  153 m  29 m (width  length  height). Therefore, the

The raw data obtained from laser scanning are smoothed by a four-order Fourier series at the pre-processing stage. The point cloud data after pre-processing are arranged into a matrix format with a size of 112  183. The knot size for the B-spline is 20. The B-spline surface model is created using 4th order B-splines in both the u and v directions. Fig. 12 plots the point cloud data and the surface model. Table 3 presents the modelling error and the execution time when the modelling algorithm is running on an Intel Q9400 Quad Core CPU with 4 GB Ram. The average error is defined as the mean absolute distance from the input data to the surface model. The results shown in Table 3 indicate that a rough layer generates more errors while the residuals decrease together with the reduction in particle size. This is mainly because the material used for the 1st layer is red scoria with an average size of 20 mm. Because scoria is a rough, vesicular, cinder-like material, and there exist many empty holes between particles after stacking, the scanning data have an undulant shape with many sharp corners (see Fig. 12a). Such sharp corner points are not differentiable in mathematical analysis. On the other hand, the nature of a B-spline surface fitting is to create a smooth surface. In mathematical terms, smooth means that the function has derivatives of all orders everywhere. Therefore, to approximate a rough surface stacked by large size particles using splines will create more errors. For the same reason, large residuals normally distribute close to the boundary of the stockpile and ground surface because the points on the edge are not differentiable. Fig. 13 plots the fitting error of the 3rd layers. A Gaussian filter smooths sharp edges on the stockpile’s surface and reduces the residuals. However, such an improvement will also produce a shrinkage effect on the 3D geometry that affects the volume calculation results. To further evaluate the modelling accuracy and assess the shrinkage in volume after smoothing, the authors introduced an isosceles triangular prism as a standard volume unit. Such primitive geometry offers a visual comparison of the modelling accuracy and gives an accurate known volume to estimate the modelling errors (because there is no direct formula

Table 4 Volume calculation of a triangular prism using QMC. unit: mm3.

Volume Error (%)

Real

Laser scanning

Surface model

Surface model (Gaussian)

53,753,375 n/a

53,654,801 0.18

52,564,307 2.21

52,526,812 2.28

Fig. 14. Surface modelling result of a prism scanned by laser.

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

12

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

Table 5 Quantity of reclaimed material for the 1st case study. Unit: g. Voxel

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

Layer 2

Table 6 Average modelling error using different degree B-spline functions. unit: cm.

Layer 3

Mass of the cut

Fe

0 0 0 0 0 0 4.10 18.72 40.64 61.90 0 0 0 0 0

0 0 0 0 0 0 2.29 11.91 23.66 39.56 0 0 0 0 0

Si 0 0 0 0 0 0 0.12 0.71 1.78 2.89 0 0 0 0 0

Al 0 0 0 0 0 0 0.10 0.34 0.90 1.81 0 0 0 0 0

Mass of the cut

Fe

43.86 44.06 58.28 70.74 82.25 93.52 100.73 95.70 83.21 72.31 52.31 65.19 79.02 93.34 102.56

24.99 25.50 33.10 40.46 49.11 54.56 63.08 58.35 46.55 37.97 23.14 26.13 34.34 41.13 50.12

Si 1.23 1.56 3.08 2.66 3.98 4.56 4.83 4.97 2.77 2.89 1.45 1.57 2.67 3.78 4.89

Layer

Degree 2nd

4th

6th

8th

10th

12th

1 2 3 4

6.32 9.48 8.98 7.98

4.92 7.61 5.64 5.08

4.79 7.39 5.47 4.42

4.39 6.31 5.09 4.15

4.36 6.21 4.87 3.96

4.31 6.01 4.49 3.87

Al 1.06 0.61 0.77 1.06 1.38 1.17 1.34 2.27 1.21 1.03 0.68 1.02 1.23 1.65 1.43

to calculate the volume of a real stockpile). Another reason to make this decision is because the triangular prism is the closest to the shape of a stockpile and is used to model a stockpile by a number of researchers in the literature, e.g. Pavloudakis and Agioutantis [5] and Lu and Maung [12]. It was made from an aluminium sheet with a dimension of 617.50 mm  174.10 mm  1000 mm (W  H  L) and scanned with the same angular resolution as the laboratory scale model, by the laser. To eliminate the errors caused by the edge detection algorithm, the boundary was identified carefully by the authors. The average dimensions of the prism observed by the author from the point cloud are

Table 7 Modelling errors of the full-scale stockpile. Layer

1 2 3 4

Number of points

B-spline model Maximum error (m)

Average error (cm)

Execution time (s)

11,076 11,076 11,076 11,076

2.18 2.04 0.80 0.72

4.99 6.91 5.08 4.14

6.25 6.22 6.27 6.31

635.0 mm  169.5 mm  997.0 mm. Points within this region were extracted from the measurement data for stockpile modelling using the same parameters. Additionally, a 3  3 Gaussian filter was applied to evaluate the shrinkage effect after the smoothing. Fig. 14 compares the surface model generated with and without the Gaussian filtering. It is clear that the surface created after the Gaussian filter has a good visual representation of the real prism. Again, the sharp edges of the prism are smoothed in the surface models.

Fig. 15. Reclaiming simulation. Phase 1, Stage 1: determine the landing points. Phase 1, Stage 2: stockpile and BWR after 10 cuts. The boundary of the 10th cut for the 3rd and 2nd layer is highlighted in cyan and black respectively. The blue dot points are the overlapped region of the third and second layer. Phase 2, Stage 1: the landing point of the second reclaim. Phase 2, Stage 2: the stockpile after the reclaim. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

13

Fig. 16. Geometric features are missing after the filtering due to the smoothing and shrinkage effects.

Fig. 17. Stockpile modelling results. The maximum error always occurs at the sharp edges or corners.

Fig. 18. Residuals of the surface model for the second layer. Most residuals are within [50 cm, 50 cm].

To calculate the volume using the QMC method, the prism region is voxelized into 70 cubes with a side length of 120 mm and each cube is filled with 65,536 low-discrepancy points. These points are substituted into the surface function for point-membership classification. As shown in Table 4, only a slight shrinkage in volume (0.07%) is observed. Also, the model generated

after the smoothing is more similar to the real object than that of without filtering. Therefore, the authors believe such smoothing filter will not introduce large errors in future the quality calculations. However, because the physical appearance of a stockpile is not solid and smooth, the authors would like to keep the Gaussian smoothing filter as an optional operation.

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

14

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

Fig. 19. Voxels of a stockpile based on the BWR reclaiming trajectory.

Fig. 20. Stockpile after the reclaiming. Regarding the cutting depth and the heights of layers, only the 3rd and 4th layer are recovered by the BWR. The geometric shapes of these two layers are plotted separately for display purposes.

Table 8 Quantity of reclaimed material for the 2nd case study. Unit: kg. Voxel

1 2 3 4 5 6 7 8 9 10

Layer 3

Layer 4

Mass of the cut

Fe

Si

Al

Mass of the cut

Fe

Si

Al

0 0 0 67.55 189.04 220.10 244.75 244.78 252.23 245.21

0 0 0 40.79 115.86 130.65 148.44 152.67 158.40 148.74

0 0 0 3.48 9.97 11.86 12.36 13.33 12.40 13.93

0 0 0 0.32 0.52 0.62 0.88 1.03 0.73 0.88

852.59 1319.41 1713.36 1811.63 1773.62 1768.30 1755.98 1737.98 1692.67 1631.36

531.43 813.74 1039.04 1082.50 1055.13 1071.85 1052.79 1055.11 999.57 979.30

39.85 64.68 86.52 97.54 99.52 78.02 76.89 86.52 100.35 90.33

2.37 6.85 10.08 9.74 4.42 5.90 7.28 2.87 9.29 5.06

The stockpile is assumed to be located on the left side of the BWR and is reclaimed by the BWR in two phases. The first phase contained 10 cuts and the second contained 5 cuts. The advancement between cuts is set to be 9 mm (0.9  (Rv  Rk)). After the first reclaim is accomplished, the BWR reverses and the BW is lowered to make the second one. To locate the potential landing points for the initial cut, the thresholds are applied, based on the BWR parameters and then a reference grid with 20 mm spacing is superimposed onto the surface model. The Gaussian smoothing filter is also applied before the B-spline interpolation. 83 candidature points are obtained after the thresholding and 38 points are validated to be effective landing points. The authors select the point with the largest cutting depth from among these validated points as the initial point for the first reclaim. Materials in the second layer are also scoped by the buckets after the 6th cut. The updated surface model is input into the auto-landing program to determine

a new landing point for the second reclaim. The states of the stockpile during the reclaiming are illustrated by a number of key frames in Fig. 15. The side length of the cube used for the QMC calculation is configured to be 30 mm. Due to the fact that the size of each cut is not identical, the side length changes slightly with the cut. The low-discrepancy sequence filled into each voxel contains 32,768 points. Based on the quality information and the voxel size, the authors divide the stockpile into 100 quality segments. The final estimation results are shown in Table 5. 4.2. Full-scale model The measurement matrix P created from the full scale dataset for B-spline fitting contains 71 rows and 156 columns. Because the geometric shape of the full size stockpile is more complex

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

but the resolution is relatively lower than the laboratory scale ones, the orders of the basic spline functions and the knot size are increased in the model. Such an increase will improve the flexibility of the B-spline function and give more freedom, allowing the spline to be bent to fit the data. With the full-scale data, the orders of basic functions in both directions are adjusted to eight and the sizes of knots are also increased to 50. These values are determined experimentally. As shown in Table 6, the improvement of the average distance error is not obvious if the degree of the basic B-spline function is greater than 8. The maximum knot size is determined by the degree of the function and the number of the row of the point data (the smaller one between the number of rows and columns of a matrix). In this case, the maximum knot size is 64 (71  8 + 1). When the size of the knot reaches 50, the maximum average error of those four layers is less than 7 cm (see Table 7), which is an acceptable value considering the size of a real stockpile. Also, no obvious improvements are observed when the size of the knot is increased continuously. The average modelling time for the entire stockpile is decreased to 6.31 s since no boundary detection is needed at the pre-processing stage. The height measurements for the ground are all zero in the point data. For the same reason, large modelling errors are mainly located around the stockpile boundaries or sharp corners. The low sampling solution may even exacerbate this situation. The maximum fitting error observed is up to 2.18 m which is at a corner of the helical-shaped cutting surface of the first layer. Table 6 shows the maximum and average error of the model. Because the second layer contains several isolated cusp-shaped regions, the average error here is the largest among the four layers. For the modelling of the full-scale data, the Gaussian smoothing filter is not suggested by the authors for two reasons: (1) The noise characteristics of the data are not known. To protect intellectual property, the data were given to the authors without detailed information. The measuring device and the accuracy of the measurement were all not provided. (2) The relatively low sampling resolution. A 2D Gaussian filter needs at least a 3  3 window to calculate the average of a point around its neighbours. Thus, smoothing is applied to a 3m2 region on the full-scale stockpile and will result in a significant loss in geometric detail, as shown in Fig. 16. Therefore, A Gaussian filter is optional if the point-to-point resolution is low. Fig. 17 illustrates the surface model of each layer separately and the corresponding distance error. The maximum error of each layer and its position are also shown in the image. Fig. 18 plots the residuals of the second layer. Although large residuals are observed in the full-scale model, the authors are still confident with the B-spline surface modelling for two reasons. First, the residuals will reduce when the point resolution is increased. Second, large residuals mainly lie in the sharp protruding or convex corners. However, only small amounts of materials exist inside these corners. Also, because the average residuals of the model are better than 7 cm, such large residuals will not cause too many errors in the quality estimations. In the reclaiming simulation, the clear distance between the BWR’s centre and the stockpile is set to be 5 m. The maximum cutting height of the BWR model is 5.88 m and the maximum height of the stockpile is 28.94 m. The gridding resolution used to generate the landing point set Pc is 0.5 m. After thresholding, 600 points on the stockpile surface are generated and 100 points pass the collision assessment. The point with a cutting depth of 17.35 m, which is very close to the effective cutting depth of the BWR (28.94 m  5.88 m = 22.06 m), is selected to activate the reclaiming operations. As the initial position of the BWR and the slewing range of the BW are obtained, the stockpile is voxelized in terms of the reclaiming trajectory (see Fig. 19). Thus, the quality and quantity of each voxel can be calculated using Eqs. (12) and (13)

15

accordingly. Fig. 20 draws the shape of the stockpile after 10 continuous cuts. The quality of recovered material in each cut is shown in Table 8. Since the real stockpile is not stacked using the chevron method, the layers are overlapped at diverse points. Such overlaps are mainly observed in the first and second layers. Although the proposed algorithm assumes the layers are topped and there are no complex intersections between adjacent layers, such regions will only create minor errors in the QMC volume calculation because the calculation is based on the surface functions. Additionally, a single cut may contain two separated 3D regions due to the shape of the folded stockpile surface. Under such circumstances, each region needs to be calculated separately. 5. Conclusion and future work In this paper, the authors present a 3D stockpile management model to assist with BWR automation and estimate the quality in blending. This is the first time that layers inside a stockpile are measured faithfully and represented in real-time. Through integrating the 3D stockpile and BWR kinematic models, the reclaiming can be controlled automatically. This is a great improvement for BWR automation because this question has not been researched in many studies. The link between these two models also allows the stockpile model to be updated in real time, according to the cutting trajectories of the BWR. Consequently, the geometric shape of the cut is controllable. Effectively, the quality of each single cut achieved by the BWR is now predictable. The chemical composition of each layer is assumed to follow normal distributions along the stacking direction. Each cut accomplished by the BWR is considered as a voxel for quality estimations. The quality of a voxel is estimated based on both its geometric information, that is, the section of each layer inside the voxel, and its variation information, that is, the normal distribution along the stacking direction. Therefore, it is possible to predict the tonnage and quality grade with a greater degree of accuracy in blending than heretofore. Based on the estimation results, it is possible to develop adaptive BWR reclaiming pattern algorithms which allow for a more flexible production, transportation and stacking schedule in iron ore handling. This 3D stockpile management model will definitely increase operational efficiency and quality control, which will benefit the entire industry. In the future, the authors expect to optimise the reclaiming pattern of the BWR so that is able to adjust the quality and approach the target grade in an effective and efficient manner, according to the quality estimation results. Additionally, if the real quality data is made available, the authors would like to develop an algorithm to simulate the quality distribution within the stockpile and use it as a database for quality prediction. Acknowledgements The authors would like to thank the Australian Research Council and Industry Partner, MatrixGroup (CMS) Pty Ltd. for their financial and in-kind support received through the ARC Linkage LP0989780 Grant. References [1] P.M. Gy, A new theory of bed-blending derived from the theory of sampling – development and full-scale experimental check, Int. J. Miner. Process. 8 (1981) 201–238. [2] J.E. Everett, Iron ore handling procedures enhance export quality, Interfaces (1996) 82–94. [3] J. Lyu, A. Gunasekaran, C.Y. Chen, C. Kao, A goal programming model for the coal blending problem, Comput. Ind. Eng. 28 (1995) 861–868. [4] G.K. Robinson, K.A. Ross, Blending in the ends of chevron stockpiles, Bulk Solids Hand. 11 (1991) 595–602.

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002

16

S. Zhao et al. / Advanced Engineering Informatics xxx (2015) xxx–xxx

[5] F.F. Pavloudakis, Z. Agioutantis, Simulation of bulk solids blending in longitudinal stockpiles, Int. J. Surf. Min. Reclam. Environ. 17 (2003) 98–112. [6] L. Tien-Fu, M. Maung Thi Rein, Optimization of reclaiming voxels for quality grade target with reclaimer minimum movement, in: 2010 11th International Conference on Control Automation Robotics & Vision (ICARCV), 2010, pp. 341–345. [7] S. Zhao, T.-F. Lu, B. Koch, A. Hurdsman, Stockpile modelling using mobile laser scanner for quality grade control in stockpile management, in: 2012 12th International Conference on Control Automation Robotics & Vision (ICARCV), 2012, pp. 811–816. [8] C. Choi, K. Lee, K. Shin, H. Ahn, Inverse kinematics of a reclaimer: redundancy and solution, in: IEEE, vol. 2883, 1997, pp. 2883–2887. [9] C. Choi, K. Lee, K. Shin, K.S. Hong, H. Ahn, Automatic landing method of a reclaimer on the stockpile, IEEE Trans. Syst., Man, Cybernet., Part C: Appl. Rev. 29 (1999) 308–314. [10] K.-H. Lee, H.-J. Bae, S.-J. Hong, Approximation of optimal moving paths of huge robot reclaimer with a 3d range finder, in: M. Gavrilova, O. Gervasi, V. Kumar, C. Tan, D. Taniar, A. Laganá, Y. Mun, H. Choo (Eds.), Computational Science and Its Applications – ICCSA 2006, Springer, Berlin/Heidelberg, 2006, pp. 151–160. [11] T.-F. Lu, M.T.R. Myo, Optimization of reclaiming voxels for quality grade target with reclaimer minimum movement, in: 2010 11th International Conference on Control Automation Robotics & Vision (ICARCV), 2010, pp. 341–345. [12] T.-F. Lu, M.T.R. Myo, Optimal stockpile voxel identification based on reclaimer minimum movement for target grade, Int. J. Miner. Process. 98 (2011) 74–81. [13] T.-F. Lu, Preparation for turning a Bucket Wheel Reclaimer into a robotic arm, in: IEEE International Conference on Robotics and Biomimetics, 2008, ROBIO 2008, 2009, pp. 1710–1715. [14] T.-F. Lu, Bucket wheel reclaimer modeling as a robotic arm, in: Proceedings of the 2009 International Conference on Robotics and Biomimetics, IEEE Press, Guilin, China, 2009, pp. 263–268.

[15] S. Zhao, T.-F. Lu, B. Koch, A. Hurdsman, 3D stockpile modelling and quality calculation for continuous stockpile management, Int. J. Miner. Process. 140 (2015) 32–42, http://dx.doi.org/10.1016/j.minpro.2015.04.012. [16] QMASTOR, Accurate stockpile management for coal and mineral mining, 2012. [17] Indurad, iStacker/iStockpile:Online 3D Heap, Dome and Stockpile/Stockyard Measurement (n.d.), 2011, . [18] S. Zhao, T.-F. Lu, B. Koch, A. Hurdsman, A simulation study of sensor data fusion using UKF for bucket wheel reclaimer localization, in: 2012 IEEE International Conference on Automation Science and Engineering (CASE), 2012, pp. 1192–1197. [19] T.-F. Lu, S. Zhao, S. Xu, B. Koch, A. Hurdsman, A 3DOF system for 3 dimensional stockpile surface scanning using laser, in: 2011 6th IEEE Conference on Industrial Electronics and Applications (ICIEA), 2011, pp. 1–5. [20] S. Zhao, T.-F. Lu, B. Koch, A. Hurdsman, Dynamic modelling of 3D stockpile for life-cycle management through sparse range point clouds, Int. J. Miner. Process. 125 (2013) 61–77, http://dx.doi.org/10.1016/j.minpro.2013. 09.009. [21] C. De Boor, A Practical Guide to Splines, Applied Mathematical Sciences, vol. 27, revised ed., Springer-Verlag, Berlin, 2001. [22] L. Piegl, W. Tiller, The NURBS Book (Monographs in Visual Communication), Springer, New York, 1997. [23] E.T. Lee, Choosing nodes in parametric curve interpolation, Comput. Aided Des. 21 (1989) 363–370. [24] T. Davies, R. Martin, Low-discrepancy sequences for volume properties in solid modelling, in: Proceedings of CSG, Citeseer, 1998, pp. 13. [25] I.M. Sobol, On the distribution of points in a cube and the approximate evaluation of integrals, USSR Comput. Math. Math. Phys. 7 (1967) 86–112. [26] D. Meagher, Geometric modeling using octree encoding, Comp. Graph. Image Process. 19 (1982) 129–147.

Please cite this article in press as: S. Zhao et al., Automatic quality estimation in blending using a 3D stockpile management model, Adv. Eng. Informat. (2015), http://dx.doi.org/10.1016/j.aei.2015.07.002