Int. J. Miner. Process. 67 (2002) 101 – 122 www.elsevier.com/locate/ijminpro
Grade of fine composite mineral particles by dual-energy X-ray radiography Gianni Schena a,*, Claudio Chiaruttini a,b, Diego Dreossi c, Alessandro Olivo c, Silvia Pani c a
Division of Georesources, Faculty of Engineering, University of Trieste, 34127 Trieste, Italy b Faculty of Natural Sciences, University of Trieste, 34127 Trieste, Italy c Department of Physics, University of Trieste, 34127 Trieste, Italy
Received 19 June 2001; received in revised form 28 February 2002; accepted 9 April 2002
Abstract Standard X-ray radiographs of a composite particle made up by two mineralogical species with different linear attenuation coefficients produce 2D attenuation images that do not allow differentiating between a thin high stopping power and a thick low stopping power-mineralogical inclusion. Dual-energy (DE) imaging—i.e., a pair of radiographs taken at different incident X-ray energy—followed by digital image subtraction allows transforming the log processed energy images into two thickness images. This linear transformation permits the estimation of the particle volumetric content in the two phases and it is fast, as required for the development of an online ‘mineral liberation’ or a ‘washability’ sensor or a ‘valuable-phase-inclusion’ detector.In this paper, micro-CT was used to reconstruct fine particle slices so that the thickness of the two phases could be precisely calculated and compared with that one derived by DE. Preliminary experiments with monochromatic X-ray demonstrated the potential of the DE method to retrieve the particle grade.The tests were conducted with fine composite particles made up of two mineralogical species. Due to the nonhomogeneous composition and density variation within the two mineralogical phases segmented and also to partial volume effects in the tomographic reconstruction of particle slices, each of the two segmented phases have a Gaussian distribution of the linear attenuation coefficient l (rather that one single value) and values of the r(l) such that the density distribution right tail of the species with the lower l overlaps the left tail of the other species giving rise to a bimodal distribution of l within the particle. This material feature adds further complexity to the problem. D 2002 Elsevier Science B.V. All rights reserved. Keywords: composite particles; mineral liberation; X-rays; dual energy
*
Corresponding author. Tel.: +39-40-676-3496; fax: +39-40-676-3497. E-mail address:
[email protected] (G. Schena).
0301-7516/02/$ - see front matter D 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 7 5 1 6 ( 0 2 ) 0 0 0 3 7 - 6
102
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
1. Motivation The fraction of composite particles of a sample that have composition within a given and narrow interval is a useful piece of information in mineral and material processing and extractive metallurgy. It can serve for the estimation of the liberation of the phases that occurs in energy-intensive size reduction processes, the evaluation of the separation efficiency of single separation units and complex concentration circuits, the calculation of the kinetics of recovery in separation/concentration processes such as flotation, leaching, etc. The following methods are/can be used for particle characterization. 1.1. Heavy-liquid methods Traditionally, fractionation by heavy liquids in classes of density has been used for the characterization of systems of particles. This method is suitable only for particles composed of two phases, for which there is a linear relationship between the average density and the volumetric composition. The laboratory procedures are tedious and use hazardous and costly chemicals (i.e., heavy liquids); the limitations of these techniques are well known and such that they are infrequently used in practice. 1.2. Stereology on polished sections Microscopic observation of the particle surface is difficult and provides biased information because the fraction of phase exposed can be different from the volumetric composition. Low-dimensional measurements of grade—i.e., the ratio of the area of valuable phase to the area of the particle transect or the ratio of length of the valuable phase to length of the sampling line within the transect—can be taken on digital images acquired with scanning electronic microscope (SEM) techniques observing a polished section made up of particles classified in narrow sizes and embodied in epoxy resin. The automated processing of these digital data with specialized image analysis techniques allows the calculation of the linear or areal grade distribution of the particle sections. The low-dimension particle grade—i.e., areal or linear ratios—distribution is a heavily distorted and meaningless representation of the volumetric grade distribution. Basically, this distortion is related to the high probability of seeing one sole phase when sectioning composite particles (these are also referred to as middling or not fully liberated particles in the processing jargon). Stereological methods, namely the kernel inversion method, are used to convert low dimension into volumetric distribution of grades (King and Schneider, 1998). 1.3. X-ray computer tomography The resolution of X-ray scanning hardware has dramatically advanced in the last decade. It is now possible to produce digital radiographs—i.e., transmission images—of the mineralogical texture of fine ore particles with a micrometric scale resolution. Thus,
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
103
the complete 3D reconstruction of the inner high-density microstructures is possible by CT. The complex mineral phase textures inside particles has been described in three dimensions with 40 Am resolution by Lin and Miller (1996). They have used reconstructions from the Michigan CT facility. Recently, these authors have commissioned a prototype cone beam CT specially developed for these mineralogical applications (Lin and Miller, 2001; Miller and Lin, 2002). However, none of the methods mentioned above appear capable of characterizing a system of particles within a time span short enough for reacting in the highly dynamic plant environments. As an example, to date, micro-CT is not suitable for on line purposes due to the time needed for a 180j rotation by small angular steps and scanning of the sample and the intense computational effort that is later required for volumetric reconstruction (unless massive dedicated parallel computing is available). 1.4. Dual-energy X-ray radiography Dual energy (DE) is a simple technique (well known in medical imaging for quantification of bone mineral content) that appears viable for the determination of the volume of the two phases that make up a composite particle. A pair of radiographs is taken at two different X-ray energies. The corresponding pair of pixels of these two digital images are ‘subtracted’ for ‘projecting out’ one phase suppressing the other. The ‘weighted subtraction’ of two energy images produces a ‘thickness’ image. The integration of the thickness images gives the volume of the projected phases. In simple words, the different attenuation properties of the two phases in the particle are used to retrieve geometric information. In the next sections, the DE theory is briefly reviewed and the experimental results are presented. DE is applied to a difficult case with two basis materials having distribution of linear attenuation coefficients with a nonnegligible variance. The DE proves to be a very promising technique for the purpose of fast estimating the grade of particles in monolayers.
2. Dual-energy theory The attenuation of the X-rays through an object is given by: I log ¼ lðEÞt; I0
ð1Þ
where l is the mean linear attenuation coefficient (cm1) of the material, I0 represents the incident (upon the object) intensity of the X-ray beam at energy E, I is the intensity emerging from the object and t is the material thickness (cm). For an X-ray path with a l varying within a nonhomogeneous material, the right side of Eq. (1) is the linear integral: mpath lðtÞdt
104
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
In this work, the interest is in two-phase particles—i.e., middlings: I log ¼ lA tA þ lB tB ; I0
ð2Þ
where subscripts A and B denote the two phases with different energy-dependent linear attenuation coefficients. 2.1. Dual-energy model Radiographing the particle twice with two different energies, one can write the following system: L ¼ M t; where: lLA M ¼ l HA
ð3Þ tA lLB log IIL0L ; L ¼ ; t ¼ ; tB lHB log IIH0H
ð4Þ
where subscripts L and H denote low and high energy, respectively. The total thickness of the two phases A and B can be resolved for: t ¼ M1 L; where: M1
ð5Þ
l 1 HB ¼ detðMÞ lHA
lLB : lLA
ð6Þ
The first (second) equation of the system (5) suppresses the effect of phase B (A) and projects out the phase A (B). A necessary condition is det(M) p 0. Large condition numbers (the ratio of the largest singular value of M to the smallest) may indicate a nearly singular matrix M that propagates small errors in the measurements of the intensities into large errors in the estimation of the thickness: this situation is avoided by proper selection of the energies. 2.2. Solution regularization with noisy energy data DE methods reportedly suffer from relatively high image noise and regularization may be necessary to eliminate the high-frequency noise: t ¼ ½MT W M þ a S1 MT W L;
ð7Þ
which minimizes the functional: 2
2
NM t LNW þ a NtNS ;
ð8Þ
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
105
where a is a coefficient that determines the relative weight given to data errors (first term) and solution regularity (second term). W is a data weighting matrix and S is a solution weighting matrix. For a=0 and W=I, the direct inversion solution is found. The default choice is: W ¼ S ¼ I;
ð9Þ
and the Tikhnov and Arsenine (1976) solution is found. A suitable value for a can be set subjectively after the visual inspection of one solution. 2.2.1. Self-tuning regularization By selecting the data weighting matrix as: 1 var ðL Þ 0 L : W ¼ 1 0 var ðLH Þ
ð10Þ
NMtLN2W =AMtLATWAMtLA has a v2 distribution. With some liberty, one can rewrite Eq. (8): 2
v2 ðaÞ þ a NtNS :
ð11Þ
The expectation for Chi-square is the degrees of freedom of the problem, m and its variance is 2m, thus, a could be set with an objective criterion: h pffiffiffiffiffii a : v2 ðaÞa mF 2m ; ð12Þ and the optimum a value is found iteratively; meaning that the data errors are allowed to be as large as the standard deviation (S.D.). These regularization methods are reviewed in some details in Press et al. (1992, Chapter 18, pp. 795 –804). The variance of L is derived in the next section (Eq. (21)). Strictly, Eq. (12) is valid with values of degree of freedom higher than 10. In such a situation, the v2 probability density is symmetric and pffiffiffiffiffican be seen as a surrogate of a Gaussian density distribution with mean m and S.D. 2m. In our case, m=2 and the v2 density distribution p(v2) is heavily skewed to the right. The following criterion was preferred: a : Pðv2 ða j m ¼ 2ÞÞ ¼ 0:99;
ð13Þ 2
where P(v2(ajv)) is the v2 probability distribution ðPðv2 Þ ¼ mv0 pðv2 Þdv2 Þ. It results to: a : v2 ðaÞ ¼ P1 ð0:99 j m ¼ 2Þ ¼ 13:81:
ð14Þ
This criterion allows as much extra regularization as statistically possible. The objective criterion also allows to check if a value of a set by visual inspection is statistically reasonable.
106
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
2.3. Further properties Further properties of the DE model (Eq. (3)) are: if tA ¼ 0
then
lHB LH ¼ ; lLB LL
ð15Þ
if tB ¼ 0
then
lHA LH ¼ ; lLA LL
ð16Þ
thus, the ratio between intensities equals the linear attenuation coefficients ratios when one phase disappears. 2.4. Variance of the thickness The overall variance of the thickness can be evaluated with simple propagation of error rules as sum of two terms, the first term accounting for the variance of the energy images and the second accounting for the variance of the linear attenuation coefficients: varðtÞ ¼ varðtÞ jvarðlÞ¼0 þvarðtÞ jvarðI =I0 Þ¼0 :
ð17Þ
2.4.1. Variance with exact linear attenuation coefficients The first term of the right side of Eq. (17) accounts for the variance of the energy images. The counting statistics of the photons that reach the detector obey the Poisson model. The detector reading and integration can alter these statistics. Here it is assumed that this added noise does not corrupt the Poisson nature of the photon counting. In an extended fashion, the thickness for the phase A is: tA ¼
1 ðl LL lLB LH Þ; detðMÞ HB
ð18Þ
thus: varðtA Þ jvarðlÞ¼0 ¼
varðtA Þ jvarðlÞ¼0 ¼
AtA ALL
2
varðLL Þ þ
1 ½detðMÞ2
AtA ALH
2 varðLH Þ
fl2HB varðLL Þ þ l2LB varðLH Þg:
ð19Þ
ð20Þ
A similar relationship can be set out for the thickness of phase B. Given the properties of the Poisson density distribution: I0 N0 var log ¼ var log ¼ varðlogN0 logN Þ I N
ð21Þ
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
var log
I0 I
¼
1 1 1 1 N þ N0 ; varðN0 Þ þ 2 varðN Þ ¼ þ ¼ N N0 N N0 N N02
107
ð22Þ
where N and N0 are the number of photons reaching the detector pixel with and without the attenuating sample. N0a[0, 2ndr1] and the same goes for N, where ndr is the dynamic range of the quantization unit. Theoretically, for an unattenuated pixel N=N0 and var (log (I0/I))=(2/N0). At the pixels attenuated by the sample, the variance is higher. 2.4.2. Variance without noise in the energy data The second term of the right side of Eq. (17) accounts for the variance of the linear attenuation coefficients. Model (5) is based on the assumption that the two phases are homogeneous, have clear interfaces and behave differently with respect to changes in the energy of the incident X-ray radiation—i.e., det(M) p 0. In practice, each phase—i.e., mineralogical species—have linear attenuation coefficients with a not negligible variance. This affects the variance of the thickness via error propagation laws: varðtA ÞjvarðI=I0 Þ¼0 ¼
" #2 X AtA Alij
ij
varðlij Þ;
ð23Þ
where:
AtA AlLA
AtA AlHA
2
AtA ¼ ; 2 Al detðMÞ LB l2HB tA2
2 ¼
2 ¼
l2HB tB2
;
ð24Þ
AtA 2 l2LB tB2 ; ¼ : detðMÞ2 AlHB detðMÞ2
ð25Þ
detðMÞ2
l2LB tA2
Thus: varðtA ÞjvarðI=I0 Þ¼0 ¼
l2HB
½tA2 varðlLA Þ þ tB2 varðlLB Þ detðMÞ2 l2LB þ ½tA2 varðlHA Þ þ tB2 varðlHB Þ; detðMÞ2
ð26Þ
and varðtB ÞjvarðI=I0 Þ¼0 ¼
l2HA
½tA2 varðlLA Þ þ tB2 varðlLB Þ detðMÞ2 l2LA þ ½t 2 varðlHA Þ þ tB2 varðlHB Þ: 2 A detðMÞ
ð27Þ
108
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
Notably, the uncertainty on the thickness of one phase is related to the value of the square of the l’s of the other phase.
3. Experimental set-up 3.1. The synchrotron radiation beamline The SYnchrotron Radiation for MEdical Physics (SYRMEP) beamline in operation at ELETTRA, the Trieste synchrotron radiation facility, was made available to our group for this experiment. The radiation source results from one of the bending magnets of the storage ring. The source dimensions (i.e., the dimensions of the electron bunches circulating in the ring) are approximately equal to 1100 (horizontal) 140 (vertical) Am2. Along the beam pipe, two beryllium windows allow the transportation of the beam from the storage ring ultrahigh vacuum to the experimental room, in air. Along the beam pipe, a monolithic channel-cut silicon (1,1,1) monochromator selects an energy band width of approximately 0.1– 0.2% in the energy range 10 – 35 keV. In this energy range, the typical flux reaching the experimental area is of about 108 photons/s mm2. In the experimental room, located at approximately 23 m from the source, a high flux, monochromatic, well-collimated laminar beam is thus available. At the entrance of the experimental room, a set of tungsten slits moved by micrometric precision stepper motors allows a precise shaping of the beam cross-section. At present, the maximum available dimensions are 100 (width) 4 (height, Full Width at Half Maximum) mm2; for this experiment, the beam cross-section was reduced down to about 30 (width) 4 (height) mm2 by means of the slit system, due to the limited detector active surface. A schematic overview of beamline layout (not to scale) is shown in Fig. 1. The imaged samples were positioned on a multi-axis, micrometric precision positioning stage. Two translation stages (accuracy=1 Am) allow sample scan in the horizontal and
Fig. 1. Elettra-SYRMEP beam line layout.
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
109
vertical directions, while a rotating stage (accuracy=103) gives the possibility of performing tomographic acquisitions. All these stages are remotely controlled via custom developed software. 3.2. The detector device A high-resolution CCD camera (X-ray HyStar Camera, Photonic Science, Millham, Mountfield, Robertsbridge, East Sussex, UK) was used as detector device. This camera features 2 2 k square pixels, each one having dimensions equal to 14 Am2; the overall active surface of the CCD camera is thus equal to 28.67 28.67 mm2. The device has a 20-Am-thick Gadolinium Oxysulphide scintillator input, placed on a straight fiber optic coupler. The sensor is cooled by a Peltier effect thermoelectric cooler, and is typically operated at 20 jC (with a coolant temperature of +10 jC). Custom developed routines allow image acquisition and manipulation, while simultaneously handling the sample positioning. This is necessary for tomographic acquisitions, where the sample is rotated by a fixed angle after each image acquisition, but desirable also when planar radiology is performed, since it allows a precise sample positioning. 3.3. Synchrotron monochromatic X-rays vs. conventional (polychromatic) sources Synchrotron radiation is characterized by a very high flux on a wide range of energies, which results in the possibility of strongly monochromatizing the beam (in the case of the SYRMEP beamline, by means of Bragg refraction from perfect crystals). Moreover, due to the extremely small source size and large source-to-sample distance, synchrotron radiation is extremely well collimated, which results a higher spatial resolution (defocusing effects due to the source size are negligible). Due to these features, synchrotron radiation is an optimal tool for DE imaging: first of all, a high flux, strictly monochromatic beam is provided, and secondly very sharp images are acquired. In principle, monochromatization via perfect crystals is feasible also when conventional sources are used, but at the price of a strongly reduced radiation flux, which results in a relevant increase of the overall acquisition time. Higher fluxes can be achieved by monochromatizing the beam provided by a conventional X-ray tube by means of mosaic crystals (see Gambaccini et al., 1995), but this flux increase would be obtained at the price of a broadened energy band width, which would result in a lower accuracy in the determination of the phase thickness. However, there is a steady improvement in traditional X-ray source technology, so that compact, inexpensive, and long-lifetime X-ray microfocus tubes are now available: microfocus tubes provide much sharper images with respect to conventional ones. Furthermore, a lot of efforts are now directed towards the development of innovative X-ray sources (see for instance Krol et al., 1997), which might in the near future provide high-intensity X-ray beams (with consequent possibility of easy monochromatization) without the necessity of relying on synchrotron radiation facilities.
110
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
4. Material Fine particles (passing 1.5 mm sieve and retained by the 1.0 mm sieve) sampled from a middling rich stream of the ‘spiral devices’ of a coal-processing plant were examined by optical microscope and those exposing the most complex phase-textures were selected for analysis (micro-CT and DE scanning). This particle size is a compromise between the desire to test fine particles and the resolution of the detector available for the experiment with a pixel of side 14 Am. However, this particle size is not to be considered the lower limit of the DE and micro-CT methods because commercial detectors with resolution of 5 Am are today available. The mean value of the linear attenuation coefficients of the marlstone matrix at various energies in the range 17– 31 keV measured on the slices reconstructed by FBP were found always comprised in the narrow band delimited by the theoretical values calculated for CaCO3 and SiO2 with PHOTCOEF (Hubbell and Seltzer, 1995/2001) software (http:// www.photcoef.com). For the carbonaceous phase, the linear attenuation coefficients were above the theoretical profile of the pure carbon at density 1.5 g ml1 (and well below the mentioned CaCO3 and SiO2 narrow band). This is because the coal phase has been contaminated by the matrix during the sedimentary rock formation process. For both the marlstone matrix and the carbonaceous phase, the mean l¯(E) are monotonic decreasing functions with respect to the energy E. While interphase limits look clear in term of luminosity on the slice image, the histogram of the pixel intensity shows that the two histograms of the individual phases overlap giving rise to an overall bimodal distribution (see Fig. 3 and next sections).
5. Methods Single particles were tomographed with a dense scanning program—i.e., 720 radiographs around 180j rotation of the turntable. Experiments were conducted in the energy range 17 –31 keV—i.e., far from the K-edges of any mineralogical element in the particle. The purpose for such a precise microtomography was two-fold. . To provide low noise slice images suitable for the segmentation of the two constituting phases (basis materials A and B) with standard image processing tools and the subsequent calculation of the phase thickness profiles to be compared to those retrieved by the DE method. . To allow an accurate determination of the linear attenuation coefficients (the elements of the matrix M) of the two phases at the specified energies used in DE experiments. Indeed, with PHOTCOEF or other specialized software, only deterministic values of the linear attenuation coefficients of the pure compounds making up the phases can be estimated. 5.1. Radiographic data processing for particle tomographic reconstruction The recorded 16-bit digital radiographs were corrected for the offset detector reading— i.e., the dark image with X-ray off—and renormalized to the image with X-ray on and no
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
111
object—i.e., the flat image—and an intensity image was obtained. Symbolically, each projection was corrected as follows: I RðhÞ D ; ¼ I0 FðhÞ D
ð28Þ
where R is the raw radiograph, D is the dark image, F is the flat image, h is the position of the turntable with respect to the starting position. The flat image (that is taken before and after the approximately 2 hour long acquisition process) was corrected for temporal variation of beam intensity—i.e., beam fluctuation— from projection to projection. The unattenuated extremes of the lines of the radiograph (‘reference channels’) at the sides of the sample were used for the line-based flat correction: ðFðh; lÞ DÞ ¼ ðFð0; lÞ DÞ
meanðRðh; l*Þ Dð0; l*ÞÞ ; meanðFð0; l*Þ Dð0; l*ÞÞ
ð29Þ
where l is the image line and l* is the line limited to the reference channels. Since the detector is 2D, many slices are acquired at the same time, each slice corresponding to one pixel row. The raw sinogram of the ith slice was obtained by orderly stacking the 720 ith rows extracted from each projection image. The sinogram was then treated for removing vertical stripes due to drifts or nonlinearities in the detector response that can give rise to ring artifact in the reconstructed image. Standard algorithms for the Radon inverse (Filtered Back Projection) assume that the axes of rotation of the sample are also the axes of the sinogram. Thus, after ring artifact reduction, the corrected sinogram was centered to overlap its central column with the projection-rotating axis. Both the baricentre of each line of the sinogram and the 0j and 180j images—i.e., the first and last projections—are used for a precise automated centering. A standard filtered back projection (FBP) algorithm was used for reconstructing the slice. The slice image is a 2D map of the distribution of the linear attenuation coefficient because the intensity of the image pixels is related to the ‘physical’ linear attenuation inside the complex mineralogical texture of the particle by: l¼
lPIXEL PIXEL SIZE
ð30Þ
where the pixel size is in centimeters. 5.2. Image analysis for phase segmentation After excluding the background, the histogram of the intensity—i.e., luminance—of the pixels of the particle slice image has the typical shape of a bimodal density distribution function. The two component functions are Gaussian density functions, each centered on the mean linear attenuation coefficient of the two mineralogical phases (l– A and l– B).
112
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
Indeed, each pixel in a CT slice represents the linear attenuation properties of a specific particle volume. Often, in a mineralized ore, that volume is comprised of a number of different substances; then, the resulting CT value represents some average of their absorption properties. As an example, the phase interfaces are not sharp and from one to the other there is a progressive change in composition reflected by the value of the linear attenuation coefficient of the corresponding voxels. The linear attenuation coefficient that corresponds to the lowest point of the ‘valley’ between the two peaks is used as threshold for segmenting the mineralogical phases. Since in our case, the peak to valley ratio was rather low, more than one slice is examined to improve the statistics of the linear attenuation coefficients. Fig. 2 shows a particle slice reconstructed by FBP from radiographs taken at 31 keV, the slice background has been changed to white to emphasize the material. The associated bimodal distribution with peaks l– PIXEL31 keV, Coal=0.00075 and l– PIXEL31 keV, Matrix=0.004 is shown in Fig. 3. 5.3. Improving the estimation of the linear attenuation coefficients The contrast in a low-energy FBP reconstructed image is higher than the one in a high-energy FBP image and the low-energy image is used for phase segmentation. The two binary image masks obtained by segmentation of the lower energy FBP intensity image are used for a better calculation the linear attenuation coefficients of the two phases. The matrix mask is shown in Fig. 4. These masks of the phases are overlapped to the particle slice images taken at the high and low energies for the precise estimation of the linear attenuation coefficients (to be used in the matrix M) ‘within’ each of the single phases, separately. The procedure is illustrated in the figures that follow for the 31 keV image. Fig. 5 show the mean values of the linear attenuation coefficient within the two phases at 31 keV. The values are calculated from the slice image in a column by column fashion. Also, the graph of the mean augmented by the standard deviation is drawn. Here, the attenuation coefficients are given in terms of lPIXEL and can be converted into cm1 by using Eq. (30). Notice that the peak of the coal phase in Fig. 3 is at a value (l– PIXEL31 keV, Coal=0.00075) lower than the calculated mean within the segmented coal phase. Vice versa, the peak of the matrix phase in Fig. 3 is at a value (l– PIXEL31 keV, Matrix=0.004) higher than the mean values with the segmented matrix phase. The same effect occurs when the values corresponding to the energy level 23 keV are taken into account. 5.4. Dual-energy imaging Some radiographs of the tomographic series were replicated at two different energies and used for DE studies. Strictly, standard DE procedures require that also Dark and Flat images are taken along with the radiograph and that each image of the Dark –Flat – Radiograph set is taken twice or more in order to improve the signal to noise ratio (SNR) by ‘frame averaging’. Temporal correction for beam decay of the white image is not necessary in standard DE experiment given that the D –F – R set of images are taken in short time sequence.
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
Fig. 2. Reconstructed slice image, 31 keV (background changed).
Fig. 3. Bimodal histogram of l (31 keV).
113
114
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
Fig. 4. Matrix binary mask.
To prove the DE method, the thickness profile of the two phases calculated with the DE algorithm was compared to the reference profile calculated on the slices (with the same orientation as the DE projections) by phase segmentation and pixel counting. An example illustrating the derivation of the reference thickness profiles is given in Figs. 6 and 7 that shows the image segmented from the original FBP and the thickness profiles measured by pixel counting, respectively. Also, the thickness of the particle (sum of the two thickness profiles) is shown. These are the reference profiles against which the DE-derived profiles are evaluated. The more adherent the DE profiles to the reference profiles the more accurate the DE method. The ‘true’ areal grade of the slice is determined by image analysis on the segmented image. The DE-derived areal grade is the integral of the DE-derived thickness profile. Symbolically, Z DE Slice Grade ¼ Z
tA ðxÞdx ;
ð31Þ
ðtA ðxÞ þ tB ðxÞÞdx
where x is the direction perpendicular to the X-ray path. Indeed, the final objective of this experiment is to prove the possibility of estimating the volumetric content of one phase in
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
115
Fig. 5. Mean and mean plus standard deviation of the linear attenuation coefficients of the two segmented phases.
Fig. 6. Segmented slice image.
116
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
Fig. 7. Reference thickness profiles.
a particle. This is also referred to as particle grade and is the sum of all the areal grades each weighted with respect to the overall slice area.
6. Experiment and results 6.1. Measurements In this section, some preliminary results of tests on single particles undertaken to demonstrate the capability of the DE method for fast particle characterization are presented. Some radiographs of the tomographic series where replicated to two different energies, thus the D – F –R sequence was not taken in short sequence and the initial D image was used. Also, the F image was derived by correcting the initial flat-image for beam decay as detailed in the previous sections. Thus, no frame averaging was applied and any image filtering avoided. A further element of disturbance that was not accounted for in DE calculation is due to the plastic straw (50 Am thick circa and 5 mm diameter circa) used to contain the particles on the rotating table and to the polystyrene-foam microspheres that filled the straw and avoided particle-particle contacts and particle movements during the tomographic acquisition. It is thought that DE measurements are penalized by these disturbing (but needed for tomography) conditions. Fig. 8 shows the raw image from FBP before the removal of the straw by Region of Interest technique.
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
117
Fig. 8. Slice before straw removal by ROI method.
The matrix M is: l23 keV;Coal M ¼ l31 keV;Coal
l23 l31
0:00250 ¼ 0:00095 keV;Matrix keV;Matrix
0:00800 10000 : 14 0:00370
ð32Þ
The coefficient (10 000/14) converts the linear attenuation coefficients lPIXEL measured on the FBP slices into l’s in (cm1) (Eq. (30)) given that CCD pixel size is 14 Am and the FBP algorithm does not apply any pixel magnification/reduction, i.e., the pixel side in the FBP images is 14 Am as in the side of the CCD pixels. Fig. 9 shows the graphs of the line radiograph—i.e., (I0/I)—at 23 and 31 keV. The higher-energy (31 keV) X-ray graph is always above the low energy one. This is because the higher the incident energy, the more penetrating the X-rays and the more transparent the object. The plastic straw and the polystyrene absorption effects are visible in the raw line radiograph. The straw edges are also responsible for the two peaks at the sides of the particle. 6.2. DE-derived areal grades on slices For the same slice used in the previous section to illustrate the experimental method, the graphs in Fig. 10 show the DE profiles obtained for a flat value of the regularization
118
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
Fig. 9. Line radiographs at 23 and 31 keV.
Fig. 10. DE-derived thickness profiles (flat a value).
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
119
Fig. 11. DE-derived thickness profiles (v2 criterion).
coefficient a=0.4 and W=I—i.e., standard Tikhonov solution. The DE-derived areal grade is 0.22 and—given the experimental conditions—in accordance with the reference value 0.27 calculated by areal measurements of the segmented image shown in Fig. 6. In Fig. 11, the objective criterion (Eq. (13)) was used to regularize the DE-derived profiles. The calculated areal grade is 0.23. It should be noticed that here the solution— i.e., the regularization of the profiles—is found with an automatic statistical criterion applied to the noisy data rather than with a visual inspection criterion forcing an aesthetic appearance of the profiles. Almost all the carried out tests indicate this level of accuracy. Finally, it should be pointed out that the matrix M is material-specific and not particlespecific. Thus, in principle, using standard DE, one can process a system of particles in a monolayer and estimate the grade of each of them in real time. Indeed, exposition time for each radiographic shot is shorter than 1 s. 6.3. Discussion Some general comments can be made: when a flat value of a is used the DEderived profile of the matrix, tB, mimics closely the matrix reference profile. Conversely, the DE-derived profile of the coal phase, tA, has rather high fluctuations. This is due to several effects, among which the higher values and the higher variance of the linear attenuation coefficients of the matrix phase and the lower thickness of the coal phase (see Eq. (27)). The uncertainty on the overall thickness is thus dominated by error on the coal thickness. However, the integral operator (Eq. (31)) has
120
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
Fig. 12. Regularization coefficient a calculated with the objective criterion.
smoothing properties and the slice grade is retrieved with a sufficient precision given that there was not previous experience on DE applied to minerals. The use of an objective regularization criterion allows to tune the value of a to the local pixel variance but this does not bring much benefits to the solution. Also, with respect to the flat regularization—i.e., a subjective and constant value of a—the objective criterion seems to over-regularize some parts of the profile and under-regularize others. This is shown by the profile of the optimum a in Fig. 12. Moreover, the procedure is more computing demanding for the optimum values of a are calculated iteratively at each pixel. The derivation of the profiles requires a calculation time of the order of seconds, a few times more time demanding than the flat regularization. By and large, one cannot conclude that the extra calculation time is paid back by an improved accuracy in the estimation of the areal grade.
7. Conclusions It appears possible to estimate with a sufficient accuracy the grade of a particle by the use of dual-energy radiography and digital subtraction methods. In the many cases examined on slices, the estimation error is often less that F15%, in spite of experimental conditions not optimized for DE (as required to prove the method against a reference slice reconstructed by microtomography). These results are very promising for they pave the way to the development of an industrial real-time sensor, returning (particle by particle) the volumetric content in the two component phases of a
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
121
monolayer of composite (and unsized) particles sampled from a stream: phase liberation or washability sensors. However, further questions need to be answered before hardware development: how is the method affected by the polychromacity of an industrial source?; how can the system maintain its performance when in duty in a severe (dusty and trembling) plant environment? On the other hand, X-ray detectors with pixel size of 5 Am are available, allowing to foresee the extension of the application to finer sizes of the particle and finer inclusion. The cost of a sensor is always to be viewed in relation to the economic benefits it brings about. However, the estimated cost of the DE hardware is competitive with that of other on line sensors today in vogue that provide much less information useful for plant control and optimization—i.e., online particle size distribution sensors. The extension of the 2 energies 2 material model (Eq. 3) to 3 energies 3 materials is not straightforward for it results in a 3 3 matrix M with a column that is linearly dependent—i.e., det(M)=0. This is because in the range of energies explored in this work, X-ray absorption by the materials is dominated by Photoelectric and Compton scattering effects: lðEÞ ¼ aCðEÞ þ bPðEÞ;
ð33Þ
where the coefficients a and b are materials dependent and independent of the energy, C(E) and P(E) are energy-dependent functions modelling the two scattering effects. Thus, the lC of the third material C may be written as a linear combination of that of the remaining two basis materials A and B. Acknowledgements The financial contributions of the Italian M.U.R.S.T. (Ministry of the University and the Scientific Research) and of the C.N.R. (National Research Council) are gratefully acknowledged. The authors like to thank Carbosulcis S.P.A. for providing the samples and the SYRMEP (Synchrotron Radiation for Medical Physics) research group. References Gambaccini, M., Taibi, A., Del Guerra, A., Frontera, F., Marziani, M., 1995. Narrow energy band X-rays via mosaic crystal for mammography application. Nucl. Instrum. Methods A 365, 248 – 254. Hubbell, J.H., Seltzer, S.M., Tables of X-ray mass attenuation coefficients and mass energy-absorption coefficients from 1 keV to 20 MeV for elements Z=1 to 92 and 48 additional substances of dosimetric interest, [Online]. Available: http://www.physics.nist.gov/xaamdi [2001, May 14]. National Institute of Standards and Technology, Originally published as NISTIR 5632, National Institute of Standards and Technology, Gaithersburg, MD (1995). King, R.P., Schneider, C.L., 1998. Stereological correction of linear grade distribution for minerals. Powder Technol. 98 (1), July 15. Krol, A. et al., 1997. Laser-based microfocused X-ray source for mammography: feasibility study. Med. Phys. 24, 725 – 732. Lin, C.L., Miller, J.D., 1996. Cone beam X-ray microtomography for three-dimensional liberation analysis in the 21st century. Int. J. Miner. Process. 47 (1 – 2), 61 – 73. June.
122
G. Schena et al. / Int. J. Miner. Process. 67 (2002) 101–122
Lin, C.L., Miller, J.D., 2001. Cone Beam X-ray microtomography—a new facility for 3-D analysis of multiphase materials preprint. SME Annual Meeting, Denver, CO. Miller, J.D., Lin, C.L., 2002. Ultimate recovery in heap leaching operation as established from mineral exposure analysis by X-ray microomography Preprint. SME Annual Meeting, Phoenix, AZ. Press, W.H., Teukolsky, A.S., Vetterling, W.T., Flannery, B.P., 1992. Integral equations and inverse theory, Chapter 18 in Numerical Recipes in Fortan, 2nd edn. Cambridge Univ. Press. Tikhonov, A.N., Arsenine, V.Y., 1976. Methodes de resolution de problemes mal poses, Translated Into French from Russian, Edition MIR, Chapter II, Section 6, pp. 73 – 79.
