Determination of Heating Value of Estonian Oil Shale by Laser-Induced Breakdown Spectroscopy

The laser-induced breakdown spectroscopy (LIBS) combined with multivariate regression analysis of measured data were utilised for determination of the heating value and the chemical composition of pellets made from Estonian oil shale samples with different heating values. The study is the first where the oil shale heating value is determined on the basis of LIBS spectra. The method for selecting the optimal number of spectral lines for ordinary multivariate least squares regression model is presented. The correlation coefficient between the heating value predicted by the regression model, and that measured by calorimetric bomb, was R2 = 0 98. The standard deviation of prediction was 0.24MJ/kg. Concentrations of oil shale components predicted by the regression model were compared with those measured by ordinary methods.


Introduction
Oil shale is used for production of fuel oil and various chemical products and is also used as a fuel in Estonian power plants.New boilers of oil shale power plants as well as oil retorts require a relatively constant quality of raw materials and fuel.The basic parameter of oil shale quality is the heating value.Heating value can vary considerably within the location in a deposit and depends substantially on concretions and limestone content [1].Optimal operation of boilers requires continuous monitoring of heating value of incoming fuel.Up to now, the heating value of oil shale in laboratories of Estonian enterprises is measured using oxygen bomb calorimeters, which are the standard instruments for measuring heating values of solid and liquid combustible samples.However, this precise method is time consuming and thus not useable for online monitoring of the continuous flow of raw material.
Laser-induced breakdown spectroscopy (LIBS) as a fast diagnostic method has been widely used for characterization of coals and lignite.Overview of activities in this field can be found for example in [2] and references therein.LIBS adapted with advanced data processing methods provides real-time information required by boiler operators nowadays.In addition to elemental composition, other important properties of coal, such as ash content and ash fusion temperature or even heating value that have complicated relation to elemental composition, are determined on the basis of LIBS spectra using multivariate analysis methods [3][4][5].Contrary to numerous studies dealing with LIBS diagnostic of coals and lignite, the number of papers related to LIBS diagnostics of oil shale is rather limited.
In [6,7], laser-induced pyrolysis was used for determination of oil shale characteristics.Directing a number of laser shots with high repetition rate to the same spot of the sample under an argon purge, the changes in the spectral line intensities in LIBS spectra over the course of measurements were related to H/C ratios determined on isolated kerogens from the rock samples.Predicted kerogen H/C ratios from the LIBS measurements of whole rock samples were well correlated (R 2 = 0 99) to values determined for kerogen isolates measured by elemental analysis [6].An important advantage of the method is that the decay of intensities of H and C lines with the growth of laser shot number allowed a clear separation of H and C in the organic part of oil shale from that in the inorganic part.However, it is not realistic to apply the above-described laser-induced pyrolysis method for online monitoring of raw material on a running belt.
In [8], elemental composition of Marcellus Shale was determined with a high precision using LIBS spectra and multivariate analysis.However, for calculation of heating value on the basis of known elemental composition, additional information about the content of elements in the organic fraction is needed.
In [9], we applied LIBS diagnostics on Estonian oil shale samples for testing the moisture-caused matrix effects.The mineral and elemental composition of tested oil shale and waste rock (kerogenic limestone) samples were very similar, except the low percentage of organics in waste rock.Changes of both, content of kerogen and content of moisture, had strong effect on intensities of analytical spectral lines.A proper choice of spectral lines allowed reliable distinction of oil shale and limestone independently of moisture content.
The oil shale chemical composition and physical properties differ significantly from those of coal and lignite.The oil shale heating value is correlated with the concentrations of elements C, H, O, S, N, and Cl in the organic fraction of oil shale [10].However, these elements are also present in large quantities in inorganic fraction of oil shale.Besides, as it is shown in [10], calculation of heating value using only elemental composition did not work well for Estonian oil shale.Mendeleev's and Dulong's formulae are widely used and accepted for calculation of heating value on the basis of elemental composition of organic part of many solid fuels, but situation with Estonian oil shale is more complicated than the applications for which these formulae were intended."There is also the possibility that describing the heat of formation of a solid fuel using only elemental composition (as is inherent in all these methods) simply does not work quite as well for oil shale organics as for other organics" [10].These reasons could complicate also the oil shale diagnostics by LIBS.Furthermore, elements belonging to surrounding air constituents cause an extra interference to LIBS spectra.According to the best of our knowledge, no investigation has been done to clarify the feasibility of LIBS for online determination of oil shale heating value.
The main aim of the present study was the evaluation of the LIBS applicability for determination of the heating value of Estonian oil shale using statistical methods of data processing.As an extra task, the LIBS-based estimation of the oil shale chemical composition was done.Experiments were done in conditions close to real industrial conditions: measurements were done in ambient air and spectra were recorded from different spots at the moving/ rotating pellet surface.About 200 g of each sample was additionally grinded using mortar and pestle to particle size smaller than 60 μm (i.e., 10 times smaller than the laser spot diameter) assuring sample homogeneity for different laser shots [11,12].Samples were pelletised without any binder, whereas before pressing 5 wt% of distilled water was added to the room dry 2 Journal of Spectroscopy powders.Powders were pressed into 25 mm in diameter and approximately 6 mm thick pellets under 18 tons held for 5 min.Four pellets from each sample were pressed.Thus, altogether 40 pellets were prepared.These pellets were stored in lab conditions at room temperature 21 °C and relative humidity about 30% during a week to balance the moisture content of samples with the humidity of ambient air resulting approximately in about 2 wt% of moisture in pellets.

LIBS Measurements.
A sketch of the laboratory LIBS setup is presented in Figure 1.The frequency-doubled Nd:YAG (YG980 Quantel) laser pulse of 532 nm wavelength and of 9 ns FWHM was focused onto the pellet surface with a plano-convex lens of 40 cm focal length in air.The focal plane of the lens was set about 3 cm under the samples' surface to avoid air breakdown in the focal point.The laser spot diameter at the sample surface was about 0.6 mm.The laser pulse energy at the plane of sample was measured with the pyroelectric laser energy meter Ophir Nova II.The pulse energy was set to 75 mJ making fluence at the pellet surface about 27 J/cm 2 .The plasma emission was collected collinearly along the laser beam axis with the off-axis parabolic mirror collimator CC52 (Andor Technology) with focal length of 52 mm and delivered via the optical fibre of 50 μm diameter to the Mechelle5000 (Andor Technology) spectrometer equipped with an intensified charge-coupled device (ICCD) camera (Andor iStar DH334T-18F-03).The spectrum in 220-850 nm range with λ/Δλ = 5000 resolution was recorded.The delay time t d = 800 ns between the laser pulse and the beginning of data acquisition was used.The recording time-gate Δt was 1 μs.These timing parameters were optimal for the signal-to-background ratio for carbon (C I) line at the wavelength 247.8 nm.Samples were fastened on a rotating stage driven with stepper motor supplying a fresh sample surface to each successive laser pulse avoiding possible pyrolysis-caused effects to the recorded spectra.The repetition rate of laser pulses was 0.55 Hz and the distance between the consecutive laser spots was more than 3 mm.Measurements were carried out under ambient atmospheric conditions.A continuous airflow of about 3 m/s was provided above the sample to remove laser-produced aerosols and dust between laser pulses [13].To suppress the effect of possible sample heterogeneity, each pellet was characterised by a sum of 100 spectra.The relative standard deviation (RSD) of the shot-to-shot fluctuations of laser energy was 0.03.

Data Processing and Results
3.1.Selection of Spectral Lines.More than a hundred emission lines belonging to 12 major elements in oil shale (C, N, H, Fe, Mg, Si, Ca, K, Al, Na, Li, and Ti) were identified in recorded spectra using NIST atomic spectra database [14].Spectral lines strongly overlapping with neighbouring lines and those with known strong self-absorption such as Ca II lines at 393.4 and 396.8 nm were removed from further analysis after preliminary investigations.
The heating value of oil shale is determined mainly by its kerogen concentration.It is known that carbon and hydrogen are the main constituent elements of kerogen [15].In 220-850 nm spectral range, the only reliably detectable carbon line was the one at 247.8 nm.In [16,17], molecular compounds C 2 and CN formed in the plasma plume were additionally used for carbon detection.Due to strong overlapping by other spectral lines, in our case, the usage of these bands was infeasible.Segments of oil shale LIBS spectra containing carbon line at 247.8 nm and hydrogen line at 656.3 nm are presented in Figure 2.
Integral intensity of a spectral line was characterised by a sum of counts belonging to pixels which cover approximately the half-width of the line (Figure 2(a)).For most of the spectral lines, the number of pixels, p, covering approximately the full width at half maximum (FWHM) comprises 9 pixels of the ICCD camera, whereas for H α line, which was strongly broadened by Stark effect, 60 pixels were used, and for whole triplet of O at 777 nm, 20 pixels were used.Usually, especially when recording with ungated detectors, for getting the actual spectral line intensities, the background should be subtracted.In the present study with gated camera and long enough delay time between the laser pulse and the recording gate, the background intensity was comparatively low.Besides, because of the presence of a number of overlapped low-intensity lines nearby some analytical lines, it was difficult to reliably measure the intensity of the background.As we checked (Section 3.4), the final result was not remarkably  influenced by background subtraction and therefore we dropped this procedure.
Poor repeatability of LIBS signal is reported in many works [12,18,19].Guidelines to reduce fluctuations of measured signals are reviewed in [12].In case of oil shale, the chemical element with a constant concentration is missing and this is why the normalisation by internal standard [18] was not applicable.
In the present study, similarly to the papers [20][21][22][23][24], the role of laser pulse energy and plasma fluctuations was diminished by normalising the line intensity with the total intensity recorded between 220 and 850 nm.
Figure 3(a) based on Table 1 shows a high correlation coefficient R 2 = 0 9846 between the mass fraction of total carbon and the heating value.At the same time, the LIBS measurements gave a correlation coefficient between the carbon line intensity and heating value only R 2 = 0 664 (Figure 3(b)).
Similarly, the correlation coefficients between the line intensities and the heating values were found for all identified spectral lines.Calculated correlation coefficients between the heating value and normalised intensities of 43 preselected lines are presented in Figure 4. Lines of low absolute value of the correlation coefficient and lines having low signal-tobackground ratio were excluded from further usage.
As one can see in Figure 4, the maximum value of correlation coefficient R = 0 861 belongs to H line. Absolute values for most of spectral lines of these correlation coefficients are less than 0.6.Since Ca and Mg lines represent mostly carbonates of waste rock, the values of correlation coefficients for Ca I, Ca II, and Mg I lines are negative.Indeed, the higher percentage of waste rock leads to a lower heating value.
It is known that one reason for the low values of correlation coefficients is matrix effect [25].Low-value correlation coefficients indicate that the prediction of heating value by only a single spectral line intensity will have low accuracy.The use of multivariate analysis helps to minimise the matrix effect of complex geochemical materials and increase the accuracy of prediction [20,26].

Multivariate Linear Regression
Model.We assumed that the heating value of an oil shale sample could be expressed as a linear combination of intensities of measured spectral lines like it was done for lignite and coal [5,27]: Here, Q denotes the heating value of a sample, vector I elements I n are intensities of spectral lines, and vector X elements x n are coefficients which values must be estimated in the course of calibration.The index n has values 0 ,…, N where N is the number of spectral lines chosen for calibration.

4
Journal of Spectroscopy Coefficients of X were calculated according to the ordinary least square regression [12,28] as where A is a M × N matrix which elements I m,n are normalised intensities of spectral lines obtained in the course of calibration and M is the number of pellets used for calibration.B is a M element vector which elements are heating values determined by calorimetric bomb method (Table 1) for samples used for calibration.Superscript T denotes the transposed matrix.

Cross-Validation.
The leave-one-out (LOO) crossvalidation method [29][30][31][32] was used to find out the optimal set of spectral lines and to estimate the regression model predictivity.We used spectra of 36 pellets with nine different heating values for calibration, and for testing, we left out a subset of spectra of four pellets of the same sample having the known heating value from calorimetric measurements.This procedure was repeated leaving out one by one each of 4-pellet subsets of 10 samples.The predictive square error (PSE) was used to evaluate the predictivity of the model: Here Q b i denotes the heating value of the ith pellet measured with bomb calorimeter (Table 1) and Q p i is calculated using relationships (1) and (2), whereby the ith spectrum is left out from the calibration set.K is the number of pellets tested (K = 40).
As the first step, PSE was calculated using only the C line for calculation of the heating value.The C line was selected because the carbon concentration is directly related to the heating value (Figure 3(a)).H α line intensity correlation with the heating value is even higher than that of C line (see Figure 4), but H α line intensity has also a strong correlation with the moisture content [33].The next step was the calculation of the heating value and PSE using additionally to the C I line one of the remaining 42 preselected lines (Figure 4).This step was repeated 42 times using one by one all remaining lines.A two-line combination with the minimal value of PSE was selected and the procedure of seeking the minimal value of PSE was repeated for combinations of 3 lines, two of them being selected from the previous step.Increasing step by step the number of lines, the minimal PSE at a certain number of lines was obtained.The dependence of minimal PSE on the number of involved lines is presented in Figure 5.
As one can see in Figure 5, the initial fast decrease of the PSE value continues until N reaches the value of about 10. Intensities of the first two lines, C and H α , in the succession of spectral lines in Figure 5 have the highest value of correlation with the heating value (Figure 4).The presence of O, Mg,   Si, Na, Ca, and Al in the initial part of progression of lines with fast decrease of PSE is also expectable as they belong to compounds presented in oil shale in relatively large quantities (Table 1, composition of samples).However, the presence of nitrogen lines there is somewhat surprising as the nitrogen content in the oil shale is very low.At the same time, nitrogen is the main component of the ambient air and its presence in the initial part of sequence could be related to the air breakdown on the surface of the samples as proposed by [34].PSE in Figure 5 has its minimum value, PSE min = 0.056, at N min = 21.Fast growth of PSE up to very high values takes place when N becomes close to the number of pellets used for calibration which is a common behaviour for ordinary least square regression.The heating value predicted by the set of first 21 lines in Figure 5 is presented in Figure 6.
The value of the correlation coefficient, R 2 = 0 981, between the heating value Q p predicted by LIBS and the calorimetric heating value Q b , is close to the correlation coefficient between Q b and the total carbon content of oil shale (Figure 3(a)).
Figure 7 presents the distribution of the difference between the "true" heating value and the heating value predicted by LIBS: Achieved standard deviation value SD = 0.24 MJ/kg differs from the uncertainty 0.13 MJ/kg of the calorimetric measurements only by a factor of two.According to Figure 5, near the vicinity of the PSE min , the value of PSE depends weakly on N. Indeed, using instead of 21 lines only the first 10 lines, gives for SD = 0.33 MJ/kg and R 2 = 0 963.

3.4.
Effect of Initial Data Variation.We tested the effect of the preselection of spectral lines as well as the procedure of the determination of integral intensity of spectral lines on the final results.Both the succession of lines and the number of spectral lines N min resulting to the minimum of PSE value (Figure 5) PSE min depended on those factors.The succession changed when for determination of line intensities, the background was taken into account, or the number of pixels, p, which was considered for integral intensity determination was varied.Similarly, the succession of lines altered with the variation in the list of preselected lines.Different successions gave a little bit different values of PSE min .Figure 8 shows how the value of PSE min and N min depend on the used numbers of pixels, p.For most of lines, the value of p was varied from 3 to 15.For wide lines O at 777 nm and H at 656 nm, the number of pixels remained unchanged.
PSE min and N min as a function of p was calculated for two preselections, one of them containing 54 lines, and another, being the subset of the first one, containing 43 lines (Figure 8).For both cases, intensity of lines was calculated with and without subtraction of the background.The background for each line was determined as the lowest intensity in spectral region 50 pixels left and right from the line.
Figure 8 demonstrates that in the case of 54 preselected lines, PSE min is always below 0.1 (MJ/kg) 2 , and in the case of 43 lines, PSE min is below 0.15 (MJ/kg) 2 .Variations of PSE min and N min with p seem to have a random character.6 Journal of Spectroscopy The shape of PSE = f(N) curves was independent of choice of p and stayed similar to that in Figure 5.The initial fast decrease of the PSE took place until N ≈ 10 followed by much slower decrease until N min .When N became close to the number of pellets used for calibration, a fast increase of PSE started.N min was below 35 in the case of 54 lines, and below 25 for 43 lines.
Like in Figure 5, when the background was not subtracted, lines of C, O, Si, Mg, H, and N were always present within the first ten lines of the succession.When the background was subtracted, the H line was not always in the list of the first ten lines.
Results presented in Sections 3.3 and 3.5 are obtained under the following conditions: (i) p = 9 because this number of pixels corresponded best to the full width at half maximum of most of spectral lines.
(ii) The background was not subtracted, because the determination of real value of background was complicated.
(iii) We use the preselection of 43 lines because the calculations are less time consuming.
However, as Figure 8 shows, most of tested values for p and other tested factors give a satisfactory (PSE min < 0.1) results as well.

Prediction of Ingredients Concentrations.
The procedure of the ordinary least square regression applied for prediction of the heating value was applied also to predict the     1).The only change introduced to the method described above was the replacing in relationships ( 1) and (3) the heating value, Q, by the concentration of a component and finding by (2) the set of coefficients which corresponds to this component.For each of the oil shale constituents determined in the laboratory, the set of spectral lines giving minimum of PSE was determined and the relative standard deviation (RSD) of prediction and correlation coefficient (R 2 ) between prediction by LIBS and laboratory measured data were calculated.Values of RSD and R 2 are presented in Table 2.
The minimum RSD = 1.5% was obtained for the total carbon.The highest values of RSD were found for constituents containing Na, K, and Al.The likely reason is that these elements were characterised by resonant lines which have a strong self-absorption [35] and as a result, the intensity of these lines is not a linear function of element concentration.
Strong sulphur lines at 921.3, 922.8, and 923.8 nm were out the wavelength range of our spectrometer.However, intensities of spectral lines of other elements are correlated with the sulphur concentration determined at the laboratory.According to Figure 9, a number of lines have values of the correlation coefficient R = 0 6 -0 8.This indicates that the presence of sulphur influences spectral line intensities of other elements.By our opinion, this influence can be attributed to the matrix effect.It was possible to predict the sulphur concentration on the basis of other lines (Figure 10) with RSD S = 4.1%.

Conclusions
The study showed that LIBS diagnostics gives reliable results for oil shale heating value in laboratory conditions using homogenised pellets of samples.
The most important results are as follows: (i) The multivariate linear regression is applicable for determination of the oil shale heating value.The used method of selection of spectral lines enables to find out minimal number of spectral lines needed for determination of the heating value and concentrations of constituents with an acceptable accuracy.
(ii) The root mean square error of predicted heating value of commercial oil shale was 0.24 MJ/kg.This value is comparable with uncertainty of determination of heating value with calorimetric bomb (0.13 MJ/kg).
(iii) The method allows to determine the dominating constituents of oil shale (C, H, CaO, SiO 2 , and Fe 2 O 3 ) with the value of relative standard deviation < 5%.
(iv) Determination of sulphur concentration is possible using a proper statistical model for description of relation between LIBS spectrum and sample composition.
stage with oil shale pellet

Figure 1 :
Figure 1: Sketch of the setup of measurements.Q-switched Nd:YAG laser YG980 and spectrometer Mechelle5000 were used.

Figure 2 :Figure 3 :
Figure 2: (a) A segment of the spectrum with carbon line (sample 1).(b) A segment with Ca and H α lines.

Figure 5 :
Figure 5: Predictive square error as a function of number of used spectral lines (N).

Figure 6 :Figure 7 :
Figure 6: Heating value of oil shale Q p determined by LIBS versus calorimetric heating value Q b .Solid lines are linear fits for data points, whereby the intercept is set to zero.

Figure 8 :
Figure 8: PSE min : minimum value of PSE as a function of number of pixels accounted for determination of integral intensities; N min : value of N at minimum PSE.

Figure 9 :Figure 10 :
Figure 9: Correlation of spectral line intensity with concentration of sulphur.

Table 1 :
Composition of dry samples.
Figure 4: Correlation coefficients between normalised spectral line intensity and heating value.

Table 2 :
RSD and R 2 for ingredients of oil shale.